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ABSTRACT 

Under this grant, two numerical algorithms were developed to predict the flow of 
viscous, hypersonic, chemically reacting gases over three-dimensional bodies. 
Both algorithms take advantage of the benefits of upwind differencing, total 
variation diminishing techniques and of a finite-volume framework, but obtain 
their solution in two separate manners. The first algorithm is a zonal, time- 
marching scheme, and is generally used to obtain solutions in the subsonic 
portions of the flow field. The second algorithm is a much less expensive, 
space-marching scheme and can be used for the computation of the larger, 
supersonic portion of the flow field. Both codes compute their interface fluxes 
with a temporal Riemann solver and the resulting schemes are made fully 
implicit including the chemical source terms and boundary conditions. Strong 
coupling is used between the fluid dynamic, chemical and turbulence 
equations. These codes have been validated on numerous hypersonic test 
cases and have provided excellent comparison with existing data. 



INTRODUCTION 


The design of recently proposed space transportation systems will, for the most 
part, be based on numerical simulations. This is because existing ground- 
based test facilities such as wind tunnels and arc jet heaters are expensive to 
operate and cannot duplicate the exact flight conditions of such vehicles. It is 
therefore the objective of this research to develop a computational tool capable 
of accurately predicting the flow of hypersonic, viscous, chemically-reacting flow 
about three-dimensional vehicles. It is the goal of this research to account for all 
relevant phenomenon in a numerical simulation. The chemical influence on a 
hypersonic flow field is apparent in the reduced heating rates to the Space 
Shuttle during reentry caused by chemical dissociation that takes place in the 
shock layer and by the non-catalytic behavior of the insulation tiles. Turbulence 
also plays an important role in performance of hypersonic vehicles and has a 
particularly strong impact on engine performance. Any proposed space 
transportation systems will also encounter these effects and hence should be 
included in the numerical simulation of such flow fields. 

Numerical solutions to steady, hypersonic flow problems can be obtained with 
both time-marching and space-marching schemes. The large amounts of 
computer time required by time-marching algorithms can prohibit the 
computation of hypersonic flow fields about realistic geometries with a 
reasonable resolution. Space-marching schemes, on the other hand, can 
provide this resolution with relatively little computer time. However, each class 
of space-marching techniques has limitations on the type of flow field it is 
capable of computing. Parabolized Navier-Stokes (PNS) algorithms are limited 
to entirely supersonic flow fields with the exception of the subsonic viscous 
layer, and Viscous Shock Layer(VSL) techniques fail in the presence of both 
axial and crossflow separation. It is therefore advantageous to employ a space- 
marching scheme in the computation of the majority of the hypersonic flow field, 
and to use a time-marching scheme only in the flow field areas where the 
space-marching scheme fails. This combination of time- and space- marching 
techniques provides the required resolution about a hypersonic configuration in 
a reasonable amount of computer time. 

Upwind numerical schemes have received wide recognition for their 
unsurpassed ability to capture flow field discontinuities without any user- 
specified smoothing terms. This property is very desirable for hypersonic 
computations in which numerous complex waves structures exist. And since 
these waves are typically very strong in such a flow field, the accurate prediction 
of their location and strength is required since small errors in these quantities 
can significantly affect the predicted dynamic loading and heat transfer on a 
vehicle. 

Thermodynamic properties can be modeled several ways in a numerical 
simulation. The perfect gas model is the most simple however is not usually 



valid for hypersonic computations. The assumption of chemical equilibrium can 
be employed in a hypersonic calculation and entails coupling a thermodynamic 
routine to the fluid dynamic equations. This can be done by either minimizing 
Gibbs Free Energy at each grid point in the computation or by utilizing a 
database of thermodynamic properties. The second approach can require much 
less computer time than the first with very little loss in accuracy. Chemical 
nonequilibrium effects are accounted for by solving additional conservation 
equations for each additional species. This can be done in either a strongly 
coupled or a loosely coupled manner. In a strongly coupled approach, the 
additional species equations are solved along with the fluid dynamic equations 
allowing information to readily flow between the two sets of governing 
equations. The weakly coupled approach allows the two equation sets to be 
solved separately which can reduce the required computer time per iteration. 
However since this second approach inhibits the flow of information between 
the two equation sets larger total computer times could become necessary to 
obtain a converged result. Finally, a thermo-chemical nonequilibrium 
computation is the most complete of the gas models by accounting for a 
nonuniform distribution of energy among the different modes. This is done by 
solving additional energy equations for each additional energy mode. 

The primary task of this grant was the development and validation of an 
accurate numerical scheme for the computation of viscous, hypersonic, 
chemically-reacting flow fields. Two codes were created for this purpose. The 
first is a zonal, time-marching scheme. It is generally used to obtain the solution 
in the subsonic or separated regions of the hypersonic flow field. Extensive 
amounts of required computer time can prohibit the computation of an entire 
hypersonic flow field with a time-marching scheme of this nature. Hence, the 
second numerical scheme was developed. It is a Parabolized Navier- 
Stokes(PNS) space-marching scheme that obtains a solution in relatively little 
computer time and can be used to compute the larger supersonic portion of the 
flow field. The time-marching code has been given the name TUFF and the 
space-marching algorithm is referred to as STUFF. TUFF stands for "A Ihree- 
Dimensional, Upwind-Differenced, fjnite-Volume Flow Solver with F.ully 
Coupled Chemistry". The additional character in the space-marching code’s 
name, STUFF, stands for Space-Marching. 

NUMERICAL SCHEMES 

The numerical schemes in both TUFF and STUFF are very similar and in fact, 
employ a lot of similar coding. This similarity is beneficial when using both 
codes on a single problem since compatibility of the two solutions is assured. It 
is also useful when making code changes since a change in one code can be 
directly implemented into the other code thereby making these codes ideal for 
research. Both codes employ a finite-volume philosophy to ensure that the 
schemes are fully conservative. Further, they obtain their upwind inviscid fluxes 
by employing a temporal Riemann solver that fully accounts for the gas model. 
A Total Variation Diminishing (TVD) technique allows extension of the schemes 
to higher orders of accuracy without introducing spurious oscillations. The 
schemes employ a strong coupling between the equation sets (including the 



fluid-dynamic, species and turbulence equations) and are made fully implicit to 
eliminate the step-size restriction of explicit schemes. This is necessary since 
step-sizes in a viscous, chemically reacting calculation can be excessively 
small for an explicit scheme, and the resulting computer times prohibitively 
large. The schemes are made implicit by fully linearizing all of the fluxes and 
source terms and by employing a modified Newton iteration to eliminate any 
numerical errors that might occur. The viscous sublayer approximation is used 
in the space-marching algorithm to allow stable marching in the presence of a 
subsonic boundary layer. Lastly, a force and moment calculation are performed 
by integrating the inviscid as well as viscous forces over all viscous boundaries. 

Both of the codes employ a versatile input routine that allows a wide range of 
problems to be solved. By changing a single parameter in the input deck, a one- 
, two- or three-dimensional analysis is possible. Further, the source terms for 
one-dimensional and two-dimensional, axi-symmetric problems are also 
options within the codes. A choice of integration options allows either an 
Approximate Facto rization(AF) or an Lower-Upper(LU) factorization of the Left- 
Hand-Side to be performed. The LU scheme has the benefit of an unlimited 
time-step size and the AF scheme contains a superior implicit boundary 
condition treatment making the choice of integration option problem-dependent. 
The codes are capable of including thin-layer viscous terms in any of the non- 
space-marching directions. An option also exists in the time-marching code to 
use a constant time-step or a constant CFL number. Finally, all grid generation 
is external to the codes allowing the user to employ state-or-the-art interactive 
grid generation procedures. 

A generalized, zonal capability was incorporated in the TUFF code to allow a 
greater flexibility of grid generation and to eliminate the memory limitations of 
large problems. It uses a patched grid technique that maintains conservation 
and allows a wide variety of grid topologies. This is done by defining the 
interface boundary as the union of all the face mesh points of the adjoining 
zones and triangulating them in a manner that maintains the original grid lines. 
Any gaps that are generated by this process are added/subtracted to/from the 
original cells in order to guarantee conservation. All of the zonal patching and 
metric computations are performed by the Conservative Interface Algorithm 
(CIA) code and are then supplied to the TUFF code for computation of the flow 
field. The CIA code serves as a preprocessor to TUFF and can also be used 
with other finite-volume codes to extend their capabilities. It has been 
successfully ported to the Compressible Navier-Stokes(CNS) code by Dr. Goetz 
Klopfer of MCAT Institute. The CIA code also has a generalized input allowing 
any boundary of any grid with any orientation to be zonal and further allows for 
curved zonal boundaries. The TUFF code uses the metrics computed by the CIA 
code to compute the inviscid and viscous interface fluxes. 

Both the TUFF and STUFF codes have generalized boundary conditions that 
allow any boundary condition to be implemented on any boundary. Boundary 
conditions must be specified for the inviscid as well as the viscous fluxes. The 
boundary conditions on the inviscid flux are listed below: 



• Zonal 

• Reflective (solid wall) 

• Non reflective (extrapolation) 

• Fixed (free stream) 

• Blunt Body Singularity 

• Continuation (no change) 

• Subsonic Inflow 

• Subsonic Outflow 

• Incompressible Inflow 

• Incompressible Outflow 

• Supersonic Blowing 


Viscous boundary conditions must be specified on velocity, temperature, and 
species gradients. These boundary condition options are listed blow: 


Velocity 


Temperature 


Catalvcitv 


• Viscous Zone 

• Inviscid 

• No Slip 

• Slip (for high altitude flows) 

• Viscous Inflow (Ablation or Blowing) 

• Viscous Outflow (Bleed) 

• Adiabatic 

• Specified Temperature 

• Radiative Equilibrium 

• Time of Flight 

• Noncatalytic 

• Fully Recombined 


For the approximate factorization option of the LHS, all of these boundary 
conditions are made fully implicit by linearizing them. The LU factorization 
option in the TUFF code only requires an explicit boundary condition. 

The codes contain several gas models including non-equilibrium, equilibrium 
and perfect gas options. An incompressible option is also included in the TUFF 
code to fill out the range of models. The choice of model is made in the input 
deck by changing a single parameter. The equilibrium option currently uses 
TannehilPs curve its for air but may be changed by modifying a subroutine for 
the thermodynamics and another for the transport properties. The 
nonequilibrium option is also easily changed. Current nonequilibrium options 
exist for air, hydro-carbon/air, hydrogen/air and the Martian atmosphere. A 
choice of backward rate computations includes specific rates or a gibbs free- 
energy minimization computation for the backward rates. Recently, a two- 
temperature model was added to a version of the TUFF code however, at the 
close of this grant period was not in working condition. Further research is 
required to debug and validate this option. 



Turbulence models that were implemented in the TUFF and STUFF codes 
include: the Baldwin-Lomax algebraic model and the Jones-Launder two- 
equation model. The Chien two-equation model was an option in early versions 
of these codes but because of problems around sharp convex corners, it was 
eliminated from the set or options. Transition is accounted for by either 
specifying the transition point or allowing the code to predict it based on Re©/M. 

For the Baldwin Lomax model, the eddy viscosity is ramped up through the 
transition region, but this is not recommended for the two-equation option. The 
Zeeman compressibility correction is an option for the two-equation model and 
improves the computation of shear layers. 

RESULTS 

Numerous validation cases were performed with the TUFF and STUFF codes 
throughout the performance of this grant. A complete list of validation cases is 
included in Appendix A. Not included in this list are any problems that were 
absent of validation data including a NASP Diverter Door study, Transonic 
Sphere study and various other studies that were of interest to NASA and the 
Principal Investigator of this grant. Several cases are highlighted below. 

Cone Test Case 

The test case that was used throughout this grant to test all of the gas models 
and other enhancements to the TUFF and STUFF codes is that of a viscous, 
hypersonic, 5-degree cone. This test case proved to be very useful since it can 
be used as a test case for both codes and since extensive data exists for 
validation. The figure below shows the results of the STUFF code with the 
chemical non-equilibrium option. 



Figure 1 . Pressure contours on a hypersonic 5 degree cone at an angle of 
attack. 



Ribbed Blunt Cone Test Case 


An experiment was conducted at the Ames Ballistic Range that provided ideal 
data for hypersonic code validation. In this experiment, blunt 5-degree cones 
were fired in air down the range. Two small shock generators (or ribs) were 
included to show the real-gas effects on the shape of the internal shock 
structure. Computations were performed with a combination of the TUFF and 
STUFF codes for validation purposes. The TUFF code was used to obtain the 
solution in the nose region and the STUFF code was used on the cone portion. 
The computed aerodynamic performance was compared with that of the 
experiment. Comparisons were also possible with the experimental 
shadowgraphs. The results of the experiment agreed quite well with those of the 
computations. 



Figure 2. Atomic oxygen contours, surface stream lines and resulting shock 
on the ribbed blunt cone geometry at M=15, a=5°. 


McDonnell Douglas Generic Option Test Case 

A hypersonic test case was performed with the McDonnell Douglas Generic 
Option geometry. Since the experimental data on this geometry was not 
available, comparison is limited to that with different CFD codes. Several 
operating conditions were computed and compared including turbulent cases. 
Below is a figure showing the resulting atomic oxygen contours of the full scale 
geometry at a Mach number of 25.3 and Reynolds number of 3,300/m. This 
figure further shows an accumulation of hot, low-speed gas on the compression 
surface that would result in a degraded performance of a scramjet engine. The 
results of this test case along with other cases at other cruise conditions were in 
good agreement with other numerical data. 



Figure 3. Atomic oxygen contours on the McDonnell Douglas Generic 
Option vehicle. 

H ypersonic Diverter Door Demonstration Case 

Of particular concern to designers of single-stage-to-orbit vehicles are the high 
aerodynamic and thermal loads on the vehicle in both the assent and reentry 
phases. A method to reduce these loads on the engine during reentry was 
conceived that included a door for the purpose of diverting the high speed flow 
around the engine. A numerical study was needed to determine the feasibility of 
this design. 

The TUFF code was used to address the feasibility of this design. Because of 
the complexity of the design, a five-zone computation was necessary. The 
computations were performed at reentry conditions (M = 15, Re = 30,000/m) 
with fully turbulent flow. The figure below shows the temperature contours of the 
TUFF results. The solution indicated that the flow was unsteady and had 
periods of high energy flow entering the cavity. These results further indicated 
that improvements were needed to continue with this design. 



Figure 4. Temperature contours on the diverter door geometry showing the 
presence of very hot gasses within the cavity. 



Hypersonic Research Vehicle Test Case 


In support of the Hypersonic Research Program, the TUFF and STUFF codes 
were utilized in the conceptual design process of a hypersonic research vehicle 
(HRV). Engineering methods are traditionally used in the conceptual design 
process but can produce erroneous results in regions where simplifying 
assumptions break down. Further, these simplified methods lack the capability 
to predict any unforeseen physics associated with a particular design. CFD, on 
the other hand, can significantly improve the accuracy and detail of the results, 
but not without a penalty. Significant computer resources can be required for a 
complete CFD analysis of a hypersonic research vehicle with an integrated 
propulsion system. 

The Ames HRV design (Mach No. = 8.0) is comprised of a waverider forebody 
with an integrated hydrocarbon scramjet engine. Waverider configurations have 
received a high degree of interest for their potentially high lift-to-drag ratios and 
their flow quality at the inlet plane. These characteristics are desirable for HRV 
missions, however the performance of the propulsion-integrated configuration 
may be much lower than that of the pure waverider shape. The propulsion 
system for the HRV is a hydrocarbon scramjet with augmented preburning. 
Hydrocarbon fuels offer sufficient specific impulse performance, heat sink 
capability and offer the potential of reduced vehicle size compared with 
hydrogen-powered designs. In addition, the handling and infrastructure 
requirements for the hydrocarbon fuels have a distinct advantage compared to 
cryogenic hydrogen. Because of the slow reaction rates of a hydrocarbon/air 
mixture in a supersonic stream, a mechanism is required to provide sufficient 
fuel/air temperatures for burning within the combustor. The concept developed 
for this purpose is an augmented preburner into which a small amount of fuel is 
burned with on-board liquid oxygen and injected into the airflow, upstream of 
the main fuel injector locations, thus ensuring that main fuel combustion is 
present and uninterrupted. 

A nose-to-tail CFD analysis was conducted on the present design to determine 
the aerodynamic performance of the integrated waverider design and to access 
the feasibility of the current engine concept. Figure 5 shows the current HRV 
design along with CFD predicted pressure contours on several cross-flow 
planes. The bow shock remains attached to the waverider leading edge for the 
entire length of the vehicle. This feature is desirable since any spillage of the 
high pressure gasses onto the upper surface would reduce vehicle 
performance. The influence of the propulsion system on the pure waverider flow 
field is also shown in this figure by numerous pressure waves. These CFD 
results also show that the flow field quality at the inlet face is very good since 
the flow is uniform and little spillage is observed. 





The current scramjet engine design is shown in the Figure 6 along with water 
contours at seven axial locations. This figure clearly shows the mechanism that 
is employed in the current design. The hot preburner gasses emerge from the 
preburner injector ports and mix with the oncoming air stream but still contain a 
very hot core just before main fuel injection station. This hot core, falling just 
above the main fuel injection, serves as a "pilot light" for main fuel injectors 
causing combustion of the main fuel to instantaneously occur. The main fuel 
injectors produce a significant amount of penetration without traversing the 
entire height of the scramjet. This figure indicates that the concept of 
preburning does indeed accomplish the task of maintaining combustion at the 
main fuel injection station and that an injector can be designed to provide 
significant flow path penetration without unstarting the engine. 


The use of CFD in the conceptual design process proved to be invaluable. It 
answered critical questions concerning the basic concepts involved in the 
design of a HRV. The nose-to-tail analysis of the waverider HRV has clearly 
shown the benefits of the current design and has revealed areas for 
improvement. The analysis of the liquid oxygen-augmented preburning 
hydrocarbon scramjet indicates that the concept is viable and does indeed 
produce uninterrupted combustion of the main fuel within the scramjet engine. 


Figure 5. 


Pressure contours about the Ames Hydrocarbon fueled waverider 
HRV. 


Onera M6 Wing Test Case 


The final test case included in this report is the Onera M6 Wing test case. This 
test case demonstrated the ability of the TUFF code to accurately predict the 
subsonic flow about a three-dimensional geometry. Comparison was made 
with experiment and with other numerical results and the TUFF results were in 
better agreement with experiment than most of the other codes. Results were 
obtained with the TUFF code with the LU option in roughly 4 hours on the Cray 
C90. The Figure below shows the surface pressure along with pressure 
comparisons with experiment. 



Chord 

Figure 8. Surface pressure contours and pressure plots at various span 
locations. 
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APPENDIX A 

TUFF/STUFF 
VALIDATION CASES 

| VALIDATED AGAINST 

CASE I EXP | CCMP l ANALY 

1) Inviscid Sphere I I I X 

2) Inviscid Ramp I I I X 

3) Viscous Flat Plate I | X | X 

4) Viscous Ramp I I X | 

5) Ram lie | X I X | 

6) Inviscid Cone I | X | 

7) Viscous Cone I I X | 

8) Sphere-Cylinder I X | X | 

9) Blunt Cone with Shock I | | 

Generators I X | j 

10) Transitional Flat Plate I X | i 

11) Transitional Cone I X | X | 

12) McDonnell Douglas Generic I I I 

Option II | X l X | 

13) Zonal Inviscid Sphere I I I X 

14) Zonal Flat Plate i j X | 

15) Zonal Viscous Sphere I | X | 

16) Transonic Nozzle I | I X 


























TUFF/STUFF VALIDATION CASES(Cont.) 


CASE 

17) Lewis Mach 5 Inlet 

18) Arc-Jet Test Case 

19) Zonal , Incompressible Cylinder 

20) Laminar, Ethelene Flame 

21) OSU Waverider Forebody 

22) Blunt Waverider 

23) Mixing/Shear Layer 

24) Ames Waverider Forebocfy 

25) Ames Hydrocarbon Scram jet 

26) Ames Hypersonic Research Vehicle 

27) Turbulent Reacting Mixing Layer 

28) Turbulent Boundary Layer 

29) Burrows-Kurkov Test Case 

30) Qnera M6 Wing Test Case 

31) Ogive Cylinder Test Case 

32) 3D Inviscid Ramp Case 


VALIDATED AGAINST | 
EXP | OCMP j ANALY | 


X | X | 


X | 


I X | 

I X I 

X | x I 

X | 
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ABSTRACT 

Two new algorithms have been developed to predict 
the flow of viscous, hypersonic, chemically reacting gases 
over three-dimensional bodies. Both take advantage of the 
benefits of upwind differencing, Total Variation Diminish- 
ing (TVD) techniques and of a finite-volume framework, 
but obtain their solution in two separate manners. The 
first algorithm is a time-marching scheme, and is gener- 
ally used to obtain solutions in the subsonic portions of the 
flow field. The second algorithm is a much less expensive, 
space-marching scheme and can be used for the compu- 
tation of the larger, supersonic portion of the flow field. 
Both codes compute their interface fluxes with a new tem- 
poral Riemann solver and the resulting schemes are made 
fully implicit including the chemical source terms. Strong 
coupling is used between the fluid dynamic and chemical 
equations. These codes have been used to compute the hy- 
personic laminar flow of reacting air over a sphere-cylinder 
and over a cone at various angles of attack. Comparison of 
the results with existing data shows good agreement. 

1. INTRODUCTION 

The design of recently proposed space transportation 
systems such as the National Aerospace Plane 1 and the 
Aeroassisted Orbital Transfer Vehicle 2 will, for the most 
part, be based on numerical simulations. This is because 
existing ground-based test facilities such as wind tunnels 
and arc jet heaters cannot duplicate the exact flight condi- 
tions of such vehicles. It is therefore necessary to model all 
relevant phenomenon in an accurate numerical simulation. 
The presence of nonequilibrium effects in a hypersonic flow 
field is apparent in the reduced heating rates to the Space 
Shuttle during reentry caused by chemical dissociation that 
takes place in the shock layer and by the non-catalytic be- 
havior of the insulation tiles. 3 The proposed space trans- 
portation systems will also encounter similar nonequilib- 
rium effects and hence these effects should be included in 
the numerical simulation of such flow fields. 
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Many efforts in the past have included various degrees 
of nonequilibrium effects in numerical simulations. Chem- 
ical nonequilibrium have been included 3-13 by solving ad- 
ditional species conservation equations in either a strongly 
coupled or a weakly coupled manner. In a strongly coupled 
approach, the additional species equations are solved along 
with the fluid dynamic equations allowing information to 
readily flow between the two sets of governing equations. 
The weakly coupled approach allows the two equation sets 
to be solved separately which can reduce the required com- 
puter time per iteration. However since this second ap- 
proach inhibits the flow of information between the two 
equation sets larger total computer times could become nec- 
essary to obtain a converged result. In addition to chem- 
ical nonequilibrium effects, thermal nonequilibrium effects 
can also be present in a hypersonic flow field. Pioneering 
research efforts 14-16 have been performed to include the 
effects thermal nonequilibrium in a numerical simulation. 
And even more recently, the effect of electro-magnetic radi- 
ation emanating from a nonequilibrium shock-layer on heat 
transfer has been addressed. 17 


Numerical solutions to steady, hypersonic flow prob- 
lems have been obtained in the past using both time-march- 
ing 4-9 and space-marching 3 ’ 10-13 schemes. The large 
amounts of computer time required by time-marching al- 
gorithms can prohibit the computation of hypersonic flow 
fields about realistic geometries with any reasonable reso- 
lution. Space-marching schemes, on the other hand, can 
provide this resolution with relatively little computer time. 
However, each class of space- marching techniques has limi- 
tations on the type of flow field it is capable of computing. 
Parabolized Navier-Stokes (PNS) algorithms 10-12 are lim- 
ited to entirely supersonic flow fields with the exception of 
the subsonic viscous layer, and Viscous Shock Layer(VSL) 
techniques 3 ’ 13 fail in the presence of both axial and cross- 
flow separation. It is therefore advantageous to employ a 
space-marching scheme in the computation of the major- 
ity of the hypersonic flow field, and to use a time-marching 
scheme only in the flow field areas where the space- marching 
scheme fails. This combination of a time-marching and a 
space-marching technique could provide the required res- 
olution about a hypersonic configuration in a reasonable 
amount of computer time. 


1 



Upwind numerical schemes have recently received wide 
recognition 18-25 for their unsurpassed ability to capture 
flow field discontinuities without any user-specified smooth- 
ing terms. This property is very desirable when comput- 
ing hypersonic flow fields where numerous complex waves 
structures can exist. And since these waves are typically 
very strong in such a flow field, the accurate prediction of 
their location and strength is required since errors in these 
quantities can drastically alter the predicted dynamic load- 
ing and heat transfer on a vehicle. 

In this paper, two numerical schemes are presented 
that incorporate several desirable features for the computa- 
tion of viscous, hypersonic, chemically-reacting flow fields. 
The first is a time-marching scheme. It is generally used 
to obtain the solution in the subsonic or separated regions 
of the hypersonic flow field. Extensive amounts of required 
computer time can prohibit the computation of an entire 
hypersonic flow field with a time-marching scheme of this 
nature. Hence, the second numerical scheme was developed. 
It is a PNS space-marching scheme that obtains a solution 
in relatively little computer time and can be used to com- 
pute the larger supersonic portion of the flow field. Both 
codes employ a finite- volume philosophy to ensure that the 
schemes are fully conservative. Further, they obtain their 
upwind inviscid fluxes by employing a new temporal Rie- 
mann solver that accounts for the presence of a multicom- 
ponent mixture of gases. A Total Variation Diminishing 
(TVD) technique of the type outlined by Chakravarthy 24 
allows extension of the schemes to higher orders of accuracy 
without introducing spurious oscillations. The schemes em- 
ploy a strong coupling between the fluid dynamic and spe- 
cies conservation equations and are made fully implicit to 
eliminate the step-size restriction of explicit schemes. This 
is necessary since step-sizes in a viscous, chemically reacting 
calculation can be excessively small for an explicit scheme, 
and the resulting computer times prohibitively large. The 
schemes are made implicit by fully linearizing all of the 
fluxes and source terms and by employing a modified New- 
ton iteration of the type described by Rai 25 to eliminate 
any linearization and approximate factorization error that 
might occur. Approximate factorization is employed to 
avoid solving many enormous banded matrices. Finally, the 
Vigneron approximation 26 is used in the space-marching al- 
gorithm to allow stable marching in the presence of a sub- 
sonic viscous layer. 


2. GOVERNING EQUATIONS 


2.1 Thin-Layer Navier-Stokes Equations 


The equations governing the flow of viscous, chemically 
reacting gases can be written in integral form as: 


d_ 

dt 



V s 


n ■ F dS = 


I 


DdV 


( 1 ) 


inviscid flux of Q through the cell faces and D consists the 
chemical source terms. In this equation, Q and D are al- 
gebraic vectors whose components are scalars and F is also 
an algebraic vector whose components are physical vectors. 

The thin-layer approximation is made by neglecting all 
viscous transport terms except those normal to a viscous 
boundary. This is justified for flow fields with the moderate 
to high Reynolds numbers. The thin-layer Navier-Stokes 
equations can now be written for the generalized, six-sided 
cell volume shown in Fig. 1 by replacing the single surface 
integral in Eq. 1 with an integral over each cell face. 

+ J J (Si+J 

V 


+ j J(F j+k -F M )d(d( + j J(G i+i - 

f(Sj+ j - Si- i m + J DdV (2) 

V 

Here, £ is the generalized streamwise coordinate, r? is the 
body normal coordinate, and ( represents the meridional 
coordinate. The indicies i,j and k represent the cell loca- 
tion in the and £ coordinate directions of the computa- 
tional mesh respectively. A non-whole index, for instance 
i + f , corresponds to a cell interface. 


By writing Eq. 1 for a generalized six-sided cell, three 
new family of fluxes are defined. The vectors E , F and G 
represent the inviscid flux through the cell interfaces with 
normals in the positive £-, rj- and ^-directions respectively. 
The dependent variables and these fluxes are presented be- 
low. 
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The metric quantities, l x , l y , and l z are components of the 
cell face normal, 1, pointing in the positive ^-direction with 
length equal to the cell face area. Similarly, the other met- 
ric quantities are the components of the vectors m and n, 
which are the normal vectors of the other family of cell 
faces, and have length equal to their area. The volume 


where V is a cell volume, n dS is a vector element of sur- 
face area with outward normal n, Q is the conserved vari- 
able per unit volume, F represents both the viscous and 
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fluxes through each family of cell interfaces can now be de- 
fined: 

U = l x u + l y v + l z w (3) 

V — m x u + m y v + m z w (4) 

W = n x u + n y v + n z w (5) 


Simply stated this expression says that the s um of the in- 
dividual species densities must equal the total density. The 
equation of state for the mixture follows Dalton’s law of 
partial pressures and is written as, 


The only viscous flux remaining in the governing equations 
after the thin-layer approximation has been made is written 
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The coefficients £%, ^ 2 > • • - ,£n +7 have been defined as 


£\~P - ml + m\ 4- m 2 z 


£2 — P j m\ + ^m 2 4- m\ 
£3 - H ml + ml 4- \rn\ 
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£& = PzpV [m 2 x 4 m 2 y + m 2 z ] 

£t+a — £&(h t — h n ) 3 = 1,2, . . . ,n — 1 

The chemical source term, D , has non-zero components 
corresponding to the species conservation equations. It is 
written as: 


where the mixture molecular weight is determined by: 


Z-2-. 


Here, c, is the species mass fraction and is defined as p t /p - 
The expression for total enthalpy is: 

5 = A + i(« , + » 2 + w ! ) = t±Z (9) 

2 p 

This expression also provides the definition of total energy, 
e. The enthalpy of the mixture is determined by summing 
the individual contributions of each species. 


fc = £c,A, 


The dimensional enthalpies and specific heats of each spe- 
cies are determined using the tables of Reference 27. Cubic 
spline interpolation is used to extract specific values from 
these tables. These tables use the following relationships 
(tildes denote dimensional quantities). 

1 . = ?<?!,,(?) + i>; ( 11 ) 

C p , - = C 2 , (12) 

where, Ci,j and C 2 , « are the values that are tabulated as 
functions of temperature, T. The frozen specific heat of the 
mixture is given by the following expression 


p/ 8 f\ 


-I ><4. 


Cl ,C 3 #= J 


where w t is the mass production rate of species 3 due to 
any number of chemical reactions involving that species. 
This production rate is dependent up on the mixture tem- 
perature and on the species molar concentrations. The ex- 
pressions for these production rates are given in Appendix 

A. 


2.2 Transport Properties 

All transport properties are those presented in Refer- 
ence 27. The viscosity of a species, 3 , is calculated using 
the following curve fit. 

P* = O.lexp [(AJog'f + B a )\o Se f + C,] (14) 


The density of the n th species is determined by the 
following mass balance: 


where A ,, B t and C a are constants for each species. Eu- 
cken’s formula is used to compute thermal conductivity 
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Wilke’s mixing rule 28 is used to compute the mixture vis- 
cosity and thermal conductivity from those of the individ- 
ual species. This is considered to. be adequate for weakly 
ionizing gases. 
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If the binary Lewis numbers for all of the species are 
assumed to be the same, then a simple expression for the 
mixture’s diffusion coefficient T> results. 


£> = 


kCe 


(19) 


Currently, the effects of multicomponent diffusion are ne- 
glected. This assumption is adequate since the molecular 
weights of the species considered are not widely different. 


2.3 Nondimensionalization 


The nondimensionalization of all the quantities pro- 
ceeds as follows (dimensional quantities are denoted by a 
tilde): 
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where Qoq is the freestream velocity and L is the geomet- 
ric length scale which is typically set to unity. Using the 
notation of Reference 11, other nondimensional quantities 
appearing in the above equations are, 
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2.4 Chemistry Model 


The chemistry model used in the present computations 
is air consisting of six species (plus electrons) undergoing 


seven reactions. The species considered are oxygen (0 2 ), 
atomic oxygen (0), atomic nitrogen (IV), nitric oxide (1V0), 
nitric oxide ion (NO- f), nitrogen (W 2 ), and electrons (e“). 
An assumption employed in this model is that the gas pos- 
sesses a zero net local charge. This allows the conservation 
of electron mass equation to be eliminated from the set of 
governing equations. The reactions that are considered are: 
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where Mi, M 2 , and M 3 are catalytic third bodies. For the 
notation in Appendix C, this model has six species, ( n = 6), 
seven reactions (m = 7) and ten reactants (n* = 10). All 
of the needed reaction rates were obtained from Reference 
27. 


3. NUMERICAL FLUX 


In this section the inviscid and viscous numerical fluxes 
are presented. First, the concept of an average flux on a 
cell interface must be addressed. This average flux(denoted 
with a “), whether it be viscous or inviscid is written as: 

# =a hII FdCd( < 2 °) 

The numerical flux, F is presented in this section and is 
an approximation of the cell face average flux, F. The 
other inviscid numerical fluxes, E and G and the viscous 
numerical flux S are represented in a similar manner. The 
average dependent variables, Q, are defined for each cell 
volume as, 


^ A^AtjAC III Qd(dr]dZ (21) 

At, A C 

The grid size parameters A£, Atj and A£ are generally 
taken to be unity. 

3.1 First Order Inviscid Flux 

The numerical inviscid flux for a first-order accurate 
upwind algorithm is determined at a cell interface with the 
help of the Riemann problem. Figure 2 shows a Riemann 
problem set up at a cell interface. In general, the states to 
the left and right of the interface are known a priori and the 
solution to the Riemann problem consists of the strength, 
direction and velocity of any waves that emanate from the 
interface. With this information, one can determine the 
invariable state remaining at the cell interface for the time 
associated with the Riemann problem. And further, one 
can determine the flux of the dependent variables through 
the cell interface . One method of obtaining an approximate 
solution to this Riemann problem is the method outlined 
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by Roe. 31 It is this method that is currently used and is 
presented below. 

In this approach, the flux at a cell interface is evaluated 
by first determining the flux change across each wave. The 
flux changes associated with waves traveling in the posi- 
tive ^-direction are given the symbol AF + and those in the 
negative direction are represented by AF~. An eigenvalue 
analysis of the matrix commonly referred to as B, re- 
veals the speed and direction of each wave. Let the matrix 
A be a diagonal matrix consisting of these wave speeds. De- 
tails of this matrix are presented in Appendix A. The flux 
remaining at the interface for all time associated with this 
Riemann problem can then be represented by the following 
equations: 

Fj+i = \(Fj + F j+J + A F~ +i - A F+ +i ) (22) 

Let R and L denote the matrix of right and left eigen- 
vectors of the matrix B respectively, evaluated at the cell 
interface. The flux difference across the positive and nega- 
tive velocity waves can be computed as follows, 

*Ff + } = 1 (r.j+i(A + |A|) i+ iL i+1 ) (Q i+1 - Qj) 

= B + (<3 j+1 - Qj) (23) 

±F- + j = i (R i+ i(A - |A|) i+ .L i+1 ) (Q i+t - Qj) 

= B-(Q j+1 - Qj) (24) 


y/Pj + y/Pi+ 1 

_ C *js/Pj + C tj+ly/Pj+l ( N 

c »;+i = n — (30) 

y/Pj + y/Pj+i 

The only remaining variables needed at the cell interface for 
the evaluation of the matricies are i , Jj+i , and D Jj+ i . 
The definition of these additional variables can be found in 
Appendix A. Since the choice of these values that satisfy the 
Roe criteria is not unique, some approximate quadratures 
are implemented. This proceeds as follows. First a good 
initial guess at these values is made. 


- _ Cj + Cj +1 

(31) 

-y _ Ti +7i+l 
1 2 

(32) 

p _ D ‘j + D ‘j + 1 

* 2 

(33) 

Next, a set of these values that satisfy Roe’s 
determined. This is accomplished by projecting 
values on to a hyperplane defined by the solution 

criteria is 
the initial 
of Eq. 25. 

D,- c 4 Ap t 6p/x 
1 — A pSp/x 

_ 7-1 +i 

7j+ s 1 _ ApSp/x 

(34) 

(35) 


where the A operator denotes [(-)j+i — (*)j| and where, 


An entropy fix is used in the definition of the upwind flux 
to avoid expansion shocks and stream line separation prob- 
lems. For further details of this entropy fix, the reader is 
referred to Reference 23. 

The matrices R, L, and A are known functions of the 
dependent variables(see Appendix B). However, the depen- 
dent variables are not defined at the cell interfaces where 
these matrices must be evaluated. Roe 21 has developed 
a special averaging procedure to calculate the dependent 
variables on the cell interface for a perfect gas and satisfy 
the following criteria, 


6p = Ap + ^ D,Ap, - (7 - 1)(A (ph) - A p) (36) 

B- 1 

n 

x = Y,( i 2 ^p-f + ( A p)* ( 37 ) 


and finally, the Roe-averaged frozen speed of sound is de- 
termined: 


V'+T 


-i 


Ti+i - x ) A j+l - Y c ‘l+i D -i+i ( 37 ) 


Fj+ 1 - Fj = (B+ + B -)%\(Qi+i - Qj) (25) 

The superscript Roc denotes a “Roe- averaged” quantity. 
Liu and Vinokur 29 have extended this analysis to include 
considerations for the presence of a multicomponent gas. 
By satisfying the relation above, the shock capturing ca- 
pabilities of the algorithm are retained and correct wave 
speeds are assured. The Roe-averaging of the dependent 
variables described by Reference 29 proceeds as follows: 


u jy/pJ Uj+ly/Pj+l 
J+ 2 y/Pj + y/Pj+l 

(26) 

v jy /p] + Vj +1 ^/p j+1 
J+ 2 y/Pj + y/Pi+1 

(27) 

W jy/Pj+ Wj+iy/pJ+i 
,J+2_ y/Pi + y/Pl^ 

(28) 


For further details of this averaging procedure, see Refer- 
ence 29. The inviscid fluxes on the families of cell interfaces 
are obtained in a similar manner by using the appropriate 
metric quantities. 

3.2 Viscous Flux 

The second-order viscous flux at the cell interface is 
obtained by realizing that all of the viscous elements have 
the form: 

- V'i) ( 38 ) 

Note that central differencing about the corresponding cell 
interface is used for any derivatives that appear in the vis- 
cous flux. If the values of <f>j + 1 represent a simple average of 
4> in the neighboring cells then a second-order evaluation of 
the viscous term results. For example in the x- momentum 
equation, the quantity that depends on the neighboring 
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3.4 Higher Order Flux 


cells and hence needs to be averaged is (()). Note that 
the cell volume V is included in the averaging. The metric 
quantities m x ,m y and m z are known quantities on the cell 
wall and hence are not included in the averaging. 

3.3 Linearization of the Fluxes 

In order to obtain either the space-marching or time- 
iterative fully implicit numerical schemes discussed later in 
this paper, certain flux Jacobians must first be determined. 
These Jacobians arise from a simple linearization of the in- 
visdd and viscous fluxes. For an implicit scheme, fluxes are 
evaluated at a location in space or time where the depen- 
dent variables are not readily known. Therefore, in order 
to evaluate these fluxes, a linearization with respect to the 
dependent variables is performed at a near by state thereby 
allowing a good approximation to this unknown flux. The 
first-order numerical flux on the j + \ cell interface, evalu- 
ated at the most recent step (whether it be a spatial step 
or time step), is represented as: 



+ (B- - - Q? +I )] (39) 

An approximate linearization of this interface flux may be 
achieved by freezing the coefficient matrix ( B~ — R + ) at 
last step p and linearizing the remaining terms. Numeri- 
cal experiments have shown that such an approximation is 
acceptable. 18 The linearized numerical flux is then written 
as: 


A higher order invisdd numerical flux can be produced 
by adding a corrective terms to the first-order flux. And in 
order to suppress any spurious oscillations that might oc- 
cur, the correction terms must fulfill certain Total Variation 
Diminishing (TVD) criteria. Reference 24 outlines a class 
of flux-difference limiting schemes that meet this criteria. 
The higher order flux (denoted by the superscript **°) is 
written as: 



(! + <£) 
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(43) 


This notation allows one to pick from a family of schemes 
with a choice of the parameter <f>. For instance, two popular 
schemes are the fully upwind scheme ( <f> — —1) and a third- 
order scheme (<j> = |). Note that the characteristic variable 
difference is the limited quantity in the present scheme. The 
characteristic variable difference is defined as: 


Aa i+ i = L i+ i (Q j+t - Qj ) (44) 


*Z\-\ fa. 

=(J?%jAQ j+I + (B L f j+i AQj + (40) 

where, A Qj is either a spatial or temporal increment in 
the dependent variables and p is the index corresponding 
to that increment. 


The limiting operators are defined with the “minmod” 
operator as, 

(•);•+! =minmod [(•),+ ! ,/?(•);_ j] ( 45 ) 

(0j+| =minmod [(•);+$,/*(•);+§] (46) 

where the minmod operator is defined as 


a< 5, = Qf 1 - 


(41) 


minmodfa y) = sign(x) * max[0, min{|x j, y * sign(x)}] (47) 


The linearization of the viscous numerical flux, S }+ 1 , 
is accomplished by freezing the value of viscosity and lin- 
earizing the remaining terms. Since these remaining terms 
are only a function of the dependent variables in the neigh- 
boring cells, the linearization is straightforward. 


4. TIME-MARCHING SCHEME 

In this section, a time-marching scheme is presented 
that can either calculate an unsteady or time-asymptotic 
solutions to chemically reacting flow problems. 


CP . ( + \ a f ", , ( 9 S I+ 1/2 

- s ^i + (ao-rJ Q,+i + \-957~ 

=-SJ + 1 + M* +i &Qj+i + Mf +i AQj 


and /? is a compression parameter that is restricted to lie 
in the range 

3 -<t> 


A similar type of linearization for the meridional flux 
G and the streamwise numerical flux E can also be per- 
formed. However, a special linearization is reserved for the 
streamwise flux in the space-marching scheme. This will be 
addressed in that section. 


4.1 Numerical Scheme 

If the high-order upwind, numerical fluxes described 
in the previous section are used at the cell faces, and an 
average set dependent variables are defined for each cell in 
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the computational mesh, the thin-layer formulation of the 
governing equations is written as: 

_ O'}. , 


At 

r>HO bHO 




\n+l 




n+1 




AHO 


n+1 


= - Vj,*) n+1 + 8) 

Note that only a first-order formula is used for the time 
term however, higher-order formulas are easily implemented 
for time-accurate computations. The numerical source term, 
D , is evaluated using the cell averaged dependent variables. 


The matrix, Tij t appearing in the above expression rep- 
resents the Jacobian of the chemical source terms that in- 
cludes derivatives with respect to both temperature and 
species concentration. Note also that A Q is defined as the 
iterative change of the cell averaged dependent variables, 
(Q p+1 — Q p ). For the first iteration, the quantities at p are 
taken from those at the time level n. Ideally, all possible 
errors are eliminated as the residual of this iteration pro- 
cess is driven to zero. However, in practice, convergence is 
defined short of this with minimal loss in accuracy. This is 
done to reduce the number of Newton iteration steps and 
hence the expense of the calculation. 


For a fully implicit scheme, the fluxes and source term 
are evaluated at the n+1 time step. This is accomplished by 
linearizing these terms in the governing equations with re- 
spect to the time- change of the dependent variables. How- 
ever, in order to limit the implicit stencil of dependent cells, 
only the first-order contribution to the inviscid fluxes is 
included in this linearization. This provides an implicit 
stencil of only three cells in each generalized coordinate 
direction. Now, to avoid the expense of inverting a large 
sparse matrix, three-dimensional approximate factorization 
is employed to break the banded matrix into three block- 
tridiagonal matrices. If the variables ■with asterisks denote 
intermediate variables, this is written in the following three 
steps. And finally, in an effort to reduce any errors that 
might occur in the scheme, such as linearization error and 
approximate factorization error, a modified Newton itera- 
tion process is employed. The entire implicit formula for 
the time-marching scheme is then written as, 
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= -(or„-o ^ 

+ (f H £ i Jt - S + (G ?? k+ 1 - G?° k 


(l-T p . k )AQ- : . k + 
M l 


At 


+ < B" - 


Re 


Vij,* 

p 


f,R * fi V A o- 
B -Tu) 


i+5 


B h - 


m r 

Re 


n+i 

p 




- [ B L - 


M l \ 


Re / 






= i.* 


5. SPACE-MARCHING SCHEME 


Space-marching schemes have been used extensively in 
the past to compute a variety of steady, supersonic flow 
fields. The numerical procedure differs from that described 
in the last section since the solution is marched in space 
rather than time. It is therefore much more efficient since 
the governing equations are only integrated once for any 
cell in the computational mesh. Further, for any cell in 
this mesh, only two (rather than three) block-tridiagonal 
inversions need be made. For these two reasons, it is easy 
to realize the benefits of having such a numerical scheme. 


The space-marching solution procedure begins by first 
providing the solver with two slabs of cells (know as the 
initial starting planes), along with all of the dependent 
variables in these cells. There are various techniques for 
obtaining this information such as the use of a either time- 
dependent code(similar to the one described above) or a 
“conical step back” procedure, 30 or even by simply spec- 
ifying freestream quantities for the entire slab. After the 
starting planes are initialized, the solution can be marched 
in the streamwise direction (^-direction) if the dependent 
variables meet a few criteria discussed in the next section. 

5.1 Parabolizing Approximation 


The mathematical nature of the governing equations 
prevents stable space-marching of the governing equations 
through subsonic or reversed flow regions of a flow field. It 
is therefore necessary to use a time-marching algorithm to 
locally obtain the solution in and around these portions of a 
hypersonic flow field. Since the viscous layer of a hypersonic 
flow field contains subsonic velocities, a spatially marched 
solution would be ill-posed and in some cases exponentially 
growing solutions (departure solutions) would be encoun- 
tered. However, a number of different techniques have been 
developed to allow stable marching of a supersonic solution 
in the presence of a subsonic viscous layer. The technique 
used in this study is that proposed by Vigneron et al. 26 


The Vigneron technique first splits the inviscid flux 
vector on the streamwise cell interface into two parts. 


E = E* + P* 


(49) 
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where, 


E* — [pU } pU U + U>l x p, pU V 4- wl V P, pUw + U)l z p, 

(t+p)V,p 1 U,p,U,...,p^ 1 U,} T 

P' = (0,(1 - u)l.P, (1 - w)i„p, (1 - u)l.p, 0, 0, 0 0] T 


The vector, E * replaces the vector, E as the inviscid flux 
through the ^-normal cell face and the vector, P* can be 
treated as a source term which is evaluated at the near- 
est supersonic point, or it can be neglected entirely. For 
this discussion and for the results presented later, it was 
neglected. 


The steady governing equations are now written by 
first eliminating the time term and by using the new stream- 
wise flux presented above: 


marching algorithm can be written in it’s entirety. 
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For a first-order accuracy in the streamwise direction, the 
streamwise numerical fluxes are evaluated as follows, 




In order to obtain a second-order streamwise flux, a second- 
order backward difference can be employed. For the results 
shown later in this paper, only a first-order difference of the 
streamwise flux was performed. 


It has been shown with an eigenvalue analysis 11 that 
stable marching in the streamwise direction is possible if: 


U3 = min 


<t 7 M| 

11 1 + (7 - l)Af| 


(52) 


where M £ is the normal component of the frozen Mach num- 
ber through the ^-normal cell interface. The factor of safety, 
<r is typically set to 0.9 

5.2 Numerical Scheme 


Finally, with the use of approximate factorization, flux 
linearization, and a modified Newton iteration, the space- 


+ { - (^X-j 

= ~ A Qi,j,k 


The matrix A* appearing in the above algorithm is the 
Jacobian of the Vigneron streamwise numerical flux. Sim- 
ilar to the time-marching scheme, A Q is defined as the 
iterative change of the cell averaged dependent variables, 
(Q p+1 — Q p ). For the first iteration, the quantities that are 
needed at i are taken as those from the last spatial step, 
i — 1. Ideally, all possible errors are eliminated as the resid- 
ual of the iteration process is driven to zero. However, as 
in the time-marching scheme, convergence is defined short 
of this with minimal loss in accuracy. This is done to re- 
duce the number of Newton iteration steps and hence the 
expense of the calculation. 


6. DECOMPOSITION 


After either a time or spatial step has been taken, all of 
the dependent variables are known at the most recent step. 
From these variables, the internal energy of the mixture, t 
can be calculated: 

£ = £ _ I( u 2 +u 2 +w 2) ( 53 ) 

p 2 

The temperature can then be obtained from the enthalpy 
curve fits with a Newton- Raphson iteration method. This 
Newton-Raphson iteration is written as, 

r‘ +1 = r* - - (54) 

5 '(T‘) ' 

where 

n 

9(T k ) = '£c.(h.( T )-hT/M.) (55) 

1=1 

g\T *) = c, (C,„(T) - MM.) (56) 

J=1 
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where the index k corresponds to the iteration index. The 
iterations are continued until convergence is reached. This 
rarely requires more than ten iterations to reduce the resid- 
ual to machine accuracy. Once the temperature is known, 
the pressure is easily computed using Eq. 7. Values for 7 
and D t are also computed at this point. 

7. FINITE VOLUME METRICS 

Before the numerical schemes detailed in previous sec- 
tions can be employed, certain metric quantities that define 
the orientation and size of each cell in the computational 
mesh must first be evaluated. These metric quantities in- 
clude three components of a surface-area vector for each cell 
interface and a volume for each cell. The generalized, six- 
sided cell shown in Fig. 1 is defined by eight verticies and 
any cell interface can be defined by the four coincident ver- 
ticies of the neighboring cells. Reference 31 gives relatively 
simple formulas for computing these metric quantities. The 
surface-area vectors are given as: 

= 2 ( r ‘+2 ~ 

m i,j+$,k — 2 ( r *~ ? *i+ 5 ,*+£ — 

x — r >- 

n »,j.*+5 ~ 2 ( r *~ 3 .i+3t*+5 ~ 

x ~ r «+|,.H-i,k+±) 

where the vectors, r represent the position vectors of the 
cell vertices. The cell volume is then computed with the 
following formula, 

( r »+ r i— i i ,A:— i ) 

There are two favorable properties of metric quantities eval- 
uated in this manner. First, the sum of all the surface-area 
vectors associated with any cell is zero. This guarantees 
that the resulting algorithm is freestream preserving. Sec- 
ond, the sum of a region of cell volumes is equal to the vol- 
ume of that region. In other words, the cell metrics allow 
no overlaping or spacing between the cells in the computa- 
tional mesh. 

8. BOUNDARY CONDITIONS 

In the finite- volume philosophy, for any boundary in 
a computational mesh, a flux must be evaluated at that 
boundary that is consistent with any boundary conditions 
that might exist. The only types of boundaries in the com- 
putation of an entirely external, supersonic flow field are: 
1) supersonic inflow(or far field), 2) supersonic outflow, 3) 
viscous boundary, and 4) symmetry conditions. Determi- 
nation of an appropriate flux for the first two boundaries is 
trivial. For a supersonic inflow the flux is simply evaluated 
using freestream quantities. And for supersonic outflow, 
the numerical scheme provides the needed flux because of 


the nature of the flow. The other two boundary fluxes are 
a little more involved. 

The symmetry condition is handled by generating ficti- 
tious cells out side of the computational domain so that the 
cell interface flux that coincides with the symmetry plane 
can be evaluated as if it were an interior face. Dependent 
variables of cells inside of the domain are reflected about 
the symmetry plane to produce a mirrored cell on the other 
side. The static dependent variables of the reflected cell ac- 
quire the identical quantity of their original. The reflected 
cell is depicted by a superscript r . 



Pm — P* s — 1, 2, . . . , n 

The dynamic variables of the reflected cell are evaluated as 
follows. 

u r = u — 2h x (n x u + n y v + h z w) 
v r — v — 2n v (n x u + n v v + n z iv) 
w r — w — 2h z (n x u + fiyV -f- h z w) 

The unit metrics in the above expression are those of a cell 
interface that is coincident with the symmetry condition. 
In order to evaluate the generalized higher-order flux at 
the symmetry condition, two reflected cells must be gener- 
ated. This symmetry condition is easily made implicit by 
realizing that the reflected cells are related to interior cells 
through this simple transformation. 

On a viscous boundary, both an inviscid and viscous 
flux must be evaluated. The inviscid flux is determined 
in a similar manner as the symmetry condition flux. The 
viscous flux is evaluated by first realizing that all of the 
velocity components are zero at the surface. Knowledge of 
either the wall temperature or heat flux must then be sup- 
plied. Evaluation of the heat flux term in the viscous flux is 
trivial if the heat flux is known a priori, however, if the wall 
temperature is known, the heat flux term is evaluated using 
a one-sided difference formula. If a catalytic wall is present, 
the mass fractions of each of the species are known at the 
wall, and the diffusion velocity terms are evaluated using 
a one-sided difference. If however, a non-catalytic wall is 
present, the diffusion velocities vanish at the wall and are 
dropped from the viscous flux entirely. 

9. RESULTS 

Two different types of geometries were chosen to vali- 
date the present techniques. The flow field about a simple 
10 degree cone is first presented at zero angle of attack 
to demonstrate the accuracy of the current methods. So- 
lutions about this same geometry at angles of attack are 
then included to show a three-dimensional capability. And 
finally, a blunt-body flow field about a sphere- cylinder ge- 
ometry is presented to demonstrate the capability of the 
current techniques to predict such flow fields. 

In this section, results obtained with both of the cur- 
rent techniques are presented. In order to distinguish the 
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results from one another, they are referenced by the names 
of their corresponding codes. The time-marching code has 
been given the name TUFF and the space-marching algo- 
rithm is referred to as STUFF. The first name stands for ” A 
Three-Dimensioned, Upwind-Differenced, Finite- Volume 
Flow Solver with Fully Coupled Chemistry”. The addi- 
tional character in the space-marching code’s name, STUFF, 
stands for Space-Marching. The similarity of the names 
of these codes reflects the similarity in the techniques and 
hence the compatability of the techniques to compute a 
single flow field. 

9.1 10° Cone Test Cases 

The first sequence of results presented here are for 
the hypersonic, laminar flow of air in chemical nonequilib- 
rium over a 10° half-angle cone at various angles of attack. 
The flow conditions were chosen to correspond to those of 
Prabhu. 11 The altitude considered was 60.96fcm and the 
free stream velocity was 8100m/s. This altitude corre- 
sponds to an ambient temperature of 252. 6K and pressure 
of 20.35iNT/m 2 and the composition by mass was assumed to 
be 26.29% molecular oxygen and 73.71% molecular nitro- 
gen. These flow conditions lead to a frozen Mach number of 

25.4 and a Reynolds number of 127, 300/m. The cone wall 
was assumed to be noncatalytic with a constant tempera- 
ture of 1200AT. For these computations, a Lewis number of 

1.4 was assumed. The angles of attack that were considered 
are 0.0, 2.5 and 5.0 degrees. 

For the 10° cone test cases presented here, results with 
both the time-marching and space- marching algorithms are 
included. The grid used for the time-marching results is 
shown in Fig. 3. The space-marching algorithm used this 
same grid as a base grid and interpolated between march- 
ing planes for any needed plane. This base grid measures 
34 cells in the streamwise direction, 30 cells in the nor- 
mal direction and 23 cells in the meridional direction. The 
length of the cone was 3.5m. The axial cell size of this base 
grid started at 0.0002m near the nose and was increased 
to nearly 0.1m at the tail end of the cone. The normal 
size of the first cell away from the body was varied linearly 
from 3 x 10~ s to 2.3 x 10 -4 . The grid was stretched from 
the body to the outer grid radius which also varied lin- 
early from 0.01m at the nose to 0.5m at the tail. For these 
test cases, all of the cross flow planes were generated to be 
axis-normal. 

The space-marching results were initiated by specify- 
ing free-stream quantities as initial data. Then at the on- 
set of the computations, a maximum CFL number of 1 
was used to determine the step size. This step size lim- 
itation was used until the bow shock emerged from the 
finely-spaced, viscous portion of the grid. The limiting CFL 
number was then increased to 30 for the remaining spatial 
steps. Using this limitation, a total of 660 spatial steps were 
required for the space-marching code. Three Newton itera- 
tions were performed for each spatial step. For the zero an- 
gle of attack test case, the required CPU-time was only 178 
seconds on the Cray-2 computer. The time-marching re- 
sults on the other hand needed much more computer time. 


The time-marching results were initiated with the space- 
marching solution and were assumed to converge after the 
maximum residual in the continuity equation decreased 3 • 
orders of magnitude. A total of about 2100 iterations were 
required. It was found that this criteria was sufficient since 
no plotable differences in the results occurred by further 
reducing this residual. For these steady-state computa- 
tions the Newton iteration process was not employed in the 
time-marching results since time- accuracy was not required. 
The resulting CPU-time needed for the zero angle of attack 
case was 3514 seconds using the time-marching technique. 
These CPU times reiterate the benefits of a space-marching 
algorithm for the computation of hypersonic flow fields. 

For the zero angle of attack case, only a two-dimensional, 
axisymmetric result was obtained. This was done by con- 
sidering only one of the meridional plane of cells in Fig. 3, 
and by assuming that the neighboring cells possess similar 
dependent variables. Figure 4 shows the axial variation of 
the surface pressure coefficient for the zero angle of attack 
case obtained with both TUFF and STUFF. Also shown in 
this figure are results of the central differencing scheme of 
Prabhu. 11 This surface pressure coefficient is defined as: 

p Pro ~ Poo 

? " iPooQlo 

All of the results predict a high pressure leading edge effect • 
that then rapidly decays to a near-constant pressure region. 
Since the Reynolds number is relatively low for this test 
case, this leading edge effect is fairly strong. Reference 11 
predicts a pressure in the later region that is about 10 per- 
cent lower than both of the current techniques. Since the 
governing equations, fluid properties and flow conditions 
of current results are identical, this discrepancy must be 
attributed to the difference in numerical procedures. Tan- 
nehill et al. 10 also found a discrepancy in the predicted 
shock shape of a centred-differencing scheme with that of 
an upwind-differencing scheme. Therefore, the discrepancy 
in Fig. 4 can in part be explained by the strong depen- 
dence of the pressure inside a shock layer on the bow shock 
shape and by the enhanced ability of upwind schemes to 
capture the bow shock. Since Reference 11 employed a 
conical-stepback procedure, this discrepancy can further 
be explained by the difference in starting solution meth- 
ods. Conical-stepback procedures make the assumption 
that the flow is locally conical however, with the specifica- 
tion of freestream at the point of the cone, this assumption 
is avoided entirely. 

Although the differences mentioned above exist be- 
tween the present results and those of Reference 11, the 
boundary layer profiles agree quite well. Figures 5-8 show 
the boundary layer profiles of various quantities at x — 
3.5 m. The abscissa of these plots corresponds to the normal 
distance, in meters, from the cone wall to the center of the 
cells. Velocity and temperature profiles are shown in Fig- 
ures 5 and 6 respectively. Excellent agreement is observed 
between the two current results, however a small discrep- 
ancy is shown near the boundary layer edge where the grid 
is relatively coarse. The peak predicted temperatures agree 
very well in all of the computations. The mass fraction pro- 
file of atomic oxygen, 0, at x = 3.5 is depicted in Fig. 7. 
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The values for mass fraction have been normalized with 
respect to the value at the wall to allow comparison with 
Reference 11. Agreement between these profiles is again ex- 
cellent. Next, the profiles of the normalized electron density 
is shown in Fig. 8. Very good agreement is observed. Fig- 
ure 9 shows the skin- friction coefficient(defined in Ref. 11) 
as a function of axial distance for the zero angle of attack 
case. Comparison of these quantities show good agreement 
with the existing data. 

In order to validate the three-dimensional capability 
of the current techniques the solution about a 10° cone at 
angles of attack is presented. The angles chosen are 2.5 
and 5.0 degrees. It was found that any greater angle re- 
quired a refinement of the base grid shown in Fig. 3. The 
flow conditions for these computations were the same as 
for the zero angle of attack case. The solution procedure, 
including the step sizes, are also identical to the previous 
case with the exception that the entire three-dimensional 
grid shown in Fig. 3 was used instead of using only a sin- 
gle meridional plane. Space- marching solutions are shown 
for both cases and time-marching solutions are only shown 
for the 2.5 degree case because of the enormity of required 
computer time for a time-marching result. 

The pressure contours in the crossflow plane at a: = 3.5 
for both angles of attack considered are shown in Figures 
10a and 10b. The quality of upwind-differenced results 
is shown in these two figures by noting that the resulting 
shock thickness does not exceed two grid cells. These fig- 
ures also depict two well known facts about cones at angles 
of attack. First, the shock layer gets thicker with greater 
angles of attack on the leeward side with a relatively small 
change in the thickness on the windward side. And second, 
it is also apparent that the pressure gradient around the 
cone becomes greater as the angle of attack is increased. In 
Figures 11a and lib, the effect on temperature with angle 
of attack is depicted. These figures present the tempera- 
ture profiles on both symmetry planes. These figures show 
that the boundary layer along with the shock layer thicken 
on the leeward side while they thin on the windward side 
with greater angles of attack. The peak temperatures how- 
ever remain unchanged with these relatively small angles of 
incidence. 

Figures 12a and 12b demonstrate the effect of angle of 
attack on chemical quantities. These figures show the O 
mass fraction profiles at an axial location of 3.5 for both 
angles of attack considered. Since the density is greater on 
the windward side of the cone, a greater degree of reaction 
can occur resulting in larger concentrations of atomic oxy- 
gen. And since the densities on the leeward side decreases 
with increasing angles of attack a reduction in the reactivity 
takes place resulting in less amounts of dissociation. These 
trends along with the thinning and thickening effects of the 
boundary layer discussed earlier are apparent in these two 
figures. Greater angles of attack are seen to only amplify 
these phenomena. Figures 13a and 13b show the effect on 
the axial variation of the heat transfer coefficient as the an- 
gle of attack is increased. The increase in the heat transfer 
rate on the windward side and a corresponding decrease on 


the leeward side with larger angles of attack is evident in 
these figures. 

9.3 Sphere-Cylinder Test Case 

This next test case is presented for two reasons. The 
first reason is to demonstrate the ability of the current tech- 
niques to predict blunt body flow fields including shock 
standoff distance. Shock standoff distance is a good mea- 
sure of the accuracy of a numerical scheme since it is di- 
rectly influenced by the distribution of dependent variables 
throughout the shock layer. The second reason that this 
test case is presented is to show how the two current tech- 
niques can be used in tandem to compute a single hyper- 
sonic flow field. The nose region of the flow field was com- 
puted with the time-marching technique, TUFF, and the 
after body was computed with the space-marching tech- 
nique, STUFF. 

The geometry and flow condition for this computa- 
tion were chosen to allow comparison with an experiment 
by Lobb 32 in which spheres were fired at hypervelocities 
into air. The bow shock locations were defined by taking 
Schlieren photographs of the spheres in flight. These condi- 
tions were also considered by Candler 16 as a test case for a 
thermo- chemical nonequilibrium code validation. Compar- 
ison of the current techniques with both references is made 
below. The flow conditions are detailed below: 

Qco = 5280(m/s) 

Toe = 293(AT) 

Poo = 664(jV/m 2 ) 

These flow conditions correspond to a frozen Mach number 
of 15.3 and a Reynolds number of 2, 190, 000/m. 

The grid used in the sphere-cylinder test case is shown 
in Fig. 14. The sphere radius was 0.635cm and the cylin- 
der portion of the geometry measured two nose radii in 
length. This grid used 45 grid cells in each normal col- 
umn of cells of which 9 are exponentially streatched near 
the body to allow for the presence of a boundary layer. 
The time-marching result used the first 24 columns of cells 
to obtain an axisymmetric solution in the nose region. The 
space-marching code used the remaining cells as a base grid 
in the computation of the afterbody region of the flow field. 
The computation proceeded by first obtaining a converged 
time-marching solution of the nose region. Machine accu- 
racy was reached after 4700 iterations at a maximum CFL 
of 10 (see Fig. 15). This required 8650 seconds on the 
Cray-2 computer. The afterbody solution was then ob- 
tained using the nose region result for the starting planes 
and marching the solution down the body. By employing 
three Newton iterations in the space-marching algorithm 
and taking 220 steps to compute the afterbody flow field, 
only 82 seconds of CPU time were required. 

Figure 16 depicts the resulting shock standoff distance 
in pressure contours about the sphere- cylinder geometry. 
Even though thermal nonequilibrium effects were neglected, 
the shock standoff distance compares very well in the entire 
nose region with the experiment of Lobb. The variation of 
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species concentration along the stagnation streamline are 
shown in Fig. 17. Comparison of the atomic oxygen profile 
with that of Candler shows good agreement. This figure 
also shows that the molecular oxygen almost completely 
dissociates in the stagnation region. Only about ten per- 
cent of the molecular nitrogen dissociates along the stag- 
nation streamline and then recombines before the wall is 
reached. This recombination is caused by the cool wall 
temperature that drives the reactions backwards. The plot 
of temperature contours in Fig. 18 also demonstrates these 
phenomena. The initial temperature rise across the strong 
bow shock in the nose region then begins to decay as a re- 
sult of the highly endothermic dissociation reactions. The 
thin thermal boundary layer is also seen in this figure. A 
contour plot of the resulting frozen Mach number contours, 
molecular oxygen contours and molecular nitrogen contours 
are shown in Figures 19, 20 and 21 respectively. 

10. CONCLUDING REMARKS 

A set of numerical schemes have been developed to 
compute the hypersonic flow of chemically reacting and 
weakly ionizing gases. These include a time-marching 
scheme and a cost effective, space-marching scheme. A fully 
implicit, strongly coupled, approximately factored method 
is employed in both techniques. The upwind, inviscid nu- 
merical fluxes are evaluated with an approximate Riemann 
solver that allows for the presence of a multicomponent 
mixture of calorically imperfect but thermally perfect gases. 
Stable higher order fluxes are obtained by utilizing a Total 
Variation Diminishing procedure. Two test cases were com- 
puted to validate the current techniques. Comparison with 
experimental data and with two existing numerical tech- 
niques shows very good agreement. Research is currently 
underway to incorporate a simple turbulence/transition 
model into the present codes. Further research includes 
adding the capability to predict hypersonic flow fields in 
various atmospheres. 
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Fig. 1 Three-dimensional Finite- Volume Cell 



Fig. 2 Illustration of a Riemann Problem 


Fig. 3 Grid Used in all Cone Computations 
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Fig. 4 Wall Pressure Coefficients 
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Fig. 5 Velocity Profiles at x = 3.5 
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Fig. 6 Temperature Profiles at x = 3.5 



0 .2 .4 .6 .8 1.0 1.2 

MASS FRACTION, c/c^,. 


Fig. 7 O Mass Fraction Profiles at x = 3.5 
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Fig. 8 Electron Molar Density Profiles at x — 3.5 
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Fig. 9 Skin Friction Coefficients 


14 





LEESIDE 

WINDSIDE 

LEESIDE 

WINDSIDE 


STUFF 


1 ®T>© 0-0 © ©-©. Q 


1 2 3 

AXIAL DISTANCE 


Fig. 13a Heat Transfer Coefficients (a — 2.5) 
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Fig. 13b Heat Transfer Coefficients (a = 5.0) 
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Fig. 14 Grid Used in Sphere- Cylinder Test Case 
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Fig. 15 Convergence Plot for CFL = 10 
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Fig. 16 Pressure Contours with Shock Standoff Comparison 
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Fig. 21 Percent iV 2 Mass Concentration Contours 


Fig. 19 Frozen Mach Number Contours 



APPENDIX A 
Inviscid Flux Jacobian 


The inviscid flux Jacobian, A, is given below: 
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where, the metrics are evaluated at the cell face and the dependent variables are those of a neighboring cell. The Inviscid 
flux Jacobians, B and C, are obtained by substituting the cell interface metrics m and n respectively in the above matrix. 

The pressure derivatives appearing in this matrix can be written as: 

dp _ _ n , (7 ~ l)(tt 2 + v* +™ 2 ) 
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3 = 1,2,... ,71 — 1 

Two convenient terms that are used in the above expressions are 7 and D,. For a mixture of calorically imperfect but 
thermally perfect gases, these terms are written as: 
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Here, C p f is the frozen specific heat of the mixture. For a single, perfect gas, D, becomes zero and 7 becomes the ratio 
of specific heats of that gas. For convenience, the temperature derivatives that appear in the viscous and source term 
Jacobians are also given. 
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APPENDIX B 
Eigenvectors 


Consider the similarity transformation between the inviscid flux Jacobian, A, and the diagonal matrix, A, consisting 
of the eigenvalues of A. This transformation is written as: 


A = LAR, L = R ” 1 


where, 

A = diag({7 + C, U - C, U, U , . . . , U) 

For this discussion, the metrics for the ^-normal cell face are used. However, the other faces are considered by simply using 
those metrics instead. The terms appearing in the eigenvalues are then defined as: 

U = l X U + lyV 4- l Z W 

C = cfo + 1J + 

where, c is the frozen speed of sound: 

A set of left and right eigenvectors providing this similarity transformation are given below. The Left eigenvectors are: 
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And the Right eigenvectors are: 
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The vector, I = l x i + l y j + l z Jc is the unit cell normal, and the vectors m and n are two arbitrary, perpendicular unit vectors 
that are in turn perpendicular to 1. The terms U , V and W are the velocities in each of these directions and can be written 
as: 

U = l x u + l y v + l z w 
V = m x u + fhyV -f m z w 
W = n x u + h y v 4- 
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APPENDIX C 
Chemical Source Terms 

The chemical source term appearing in Eq. 1 is written as: 

D = (0,0,0, 0,0, w t ,w2,...,w n - 1 ) T 

where w, is the mass production of species s resulting from any number of chemical reactions involving that species. 
Consider a mix ture of gases with n species undergoing m simultaneous reactions. Each chemical reaction is represented 
symbolically as: 

Y. v kiAi ^ y, v'lm 

i=i i=i 

where u’ kl and i are the stoichiometric coefficients and Ai is the chemical symbol for the l tfl species, rit is the total number 
of reactants including third bodies and electrons. The mass production rate of species s is determined by the law of mass 
action. It is stated as: 

m [" rv, tv, 

w, = M t ”, - I'i'.J k ft l - h,l UM- 

1=1 L r= 1 r=l 

The molar concentrations, 7 ,., of the species are defined as: 

7r — Pr / AA r T = 1,2,... ,n 

and the molar concentrations of the catalytic third bodies are expressed in terms of their third body efficiencies as: 

n 

7r = ^(r-*),. 7 « r = n + 1 , n + 2, . . . , (n + n tb ) 

•= \ 

where ntb is the number of third bodies considered. If electrons are considered, their molar concentration is evaluated by 
assuming the molar concentration of electrons is equal to the molar concentration of ions. This is stated as, 


7n, 


, = X^«7. 


where, I t is the charge of species s and the subscript n* represents the last considered species which in this case is the 
electron. 

The forward and backward reaction rates of any reaction considered are assumed to be solely a function of temperature 
and are expressed in the modified Arrhenius form as 
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These reaction rate coefficients have been nondimensionalized with the following expressions. 
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The terms C 3< i, C 2 ,i, &3,l, &i ,i» -^ 2,1 and D 3y i are constants for a particular reaction l. 
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EXPERIMENTAL AND COMPUTATIONAL RESULTS FOR 5 DEGREE 
BLUNT CONES WITH SHOCK GENERATORS AT HIGH VELOCITY 


by 

A. W. Strawa*, G. A. Molvik**, L. A. Yates* , and C. Comelison* 
NASA Ames Research Center, Moffett Field, California 


ABSTRACT 


Experiments and computations have been performed under laminar conditions in air on 5° blunt cones at 
velocities of 5 km/s and 6 km/s and at Reynolds numbers of 10 5 and 10 6 . The experiments were conducted in the 
Ames’ ballistic range. The computations were performed using ideal-, equilibrium- and nonequilibrium-chemistry 
models for air. At the conditions of the tests, the aerodynamic coefficients are sensitive to the real-gas effects present, 
and both experimental and computational aerodynamic coefficients show real-gas and non-linear effects. The nonequi- 
librium computations show that a large amount of oxygen is dissociated in the blunt nose region of the flow and 
much of the oxygen remains dissociated over the entire length of the body, providing an insight into the source of 
the observed effects in the aerodynamic coefficients. The experimental and computational shock-shapes are in good 
agreement. 


NOMENCLATURE 


A 

C D 

C l 

Cm 

C Md 


d 

Moo 

m 

Rei 


T 

(u,v,w) 

V 


(x,y,z) 

a 

OCT 

P 

e 

7 

(Z 

Q 

q* 

p 


model base area 
drag coefficient 
lift coefficient 
pitching moment coefficient 
combined damping coefficient 
Cmj. - Cm, + Cm„ 
model base diameter 
free-stream Mach number 
model mass 

free-stream Reynolds number based on 

model length 

temperature 

model velocity with respect to body axis 

velocity 

V — vu 2 + v 2 + w 2 

model position in Cartesian system 

pitch angle 

thermal diffusivity 

yaw angle 

estimated error in aerodynamic coefficients 

ratio of specific heats 

thermal conductivity 

rate of pitch 

surface heat flux 

density 
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a 

subscripts 

resultant angle-of-attack; 
sin(cr) = s/v 1 + w 1 jV 

b 

base 

n 

nose 

3 

surface value 

oo 

free stream value 


INTRODUCTION 

Future missions of national interest require the 
design of vehicles that will fly at hypervelocities and at 
high altitudes. 1,2 New mission profiles will subject these 
future vehicles to harsh environments where not all of 
the relevant physical phenomenon are completely un- 
derstood. Since no ground based facility can completely 
simulate the flight environment, computational solutions 
will become increasingly important in vehicle design 
and in predicting vehicle performance. To this end, it is 
essential that our ground based facilities and computer 
codes be used together to understand the relevant phe- 
nomenon and to insure that they are properly modeled. 
Experiments that can better isolate the effects of cer- 
tain physical phenomena can be designed with the aid of 
computational aerothermodynamics, and experimental re- 
sults can be used to verify the computational models and 
methods. The process of code verification is discussed 
in Ref. 3. 

The objective of this research was to conduct 
an experiment which would produce non-linear and real- 
gas effects in the aerodynamic coefficients for a generic 
model, to use a state-of-the-art computer code to com- 
pute the flowfields, and to compare the results. Past ex- 
perience with the Apollo vehicle and the Space Shuttle 
(STS) indicate that real-gas effects can have a signifi- 
cant influence on pitching moments. 4-6 Computational 



methodology at that time was unable to predict the full 
extent of these effects. Computational methods today 
have matured to the point were complex hypersonic 
flowfields can be calculated. However, experimental 
data demonstrating real-gas effects on slender bodies 
which are suitable for code verification are very lim- 
ited. Welsh et al. conducted tests on 10° blunt cones at 
V = 4 .8 km/ s in nitrogen and air in a ballistic range. 7 
They attributed changes in Cm to real-gas effects. Kruse 
obtained drag data for 6° blunt cones at V = 5 2km/ s 
and Rei = 300,000 in the Ames’ ballistic range. 8 Mal- 
colm and Rakich showed the effect of nose bluntness 
on 12.5° cones at Mach 17. 9 All of these tests showed 
nonlinearity and real-gas effects in the aerodynamic co- 
efficients. Real-gas effects were observed by Strawa et 
al. on sharp cones. 10 Additionally, real-gas effects have 
been simulated by matching 7 in the Langley CFa tun- 
nel (see, for example, ref. 11). 

In the present study, a 5° half-angle cone with 
a bluntness ratio, d*. /db, of .3 was selected. Two small 
shock generators were included to show the effects of 
real-gas chemistry on the shape of internal shocks. The 
velocities were 5 and 6 km / 3 and the Mach numbers 
were 14 .5 and 17 .6 . The computations were conducted 
using a combination of time- and space- marching algo- 
rithms which included fully-coupled chemistry models. 
The strong shock produced by the blunt nose at these 
speeds generated temperatures high enough to dissociate 
a significant amount of oxygen. The Reynolds num- 
bers (based on model length) were varied from 10 6 to 
10 5 , which should result in a increase in skin friction 
by a factor of three for laminar flow conditions. Viscous 
interactions between the boundary layer and the outer 
inviscid flow can have important effects on the surface 
pressure distribution and, hence, on the drag, lift, and 
moment coefficients. 12 The moment coefficient of the 
present model should be sensitive to real-gas effects on 
the pressure and skin friction. 

THE BALLISTIC RANGE EXPERIMENTS 

The experimental data presented in this pa- 
per were obtained in the the Hypervelocity Free-Flight 
Aerodynamic Facility (HFFAF) ballistic range at Ames 
Research Center. The HFFAF ballistic range is shown 
schematically in Fig. 1. The facility consists of: a two- 
stage light-gas gun, a tank in which the sabot is sepa- 
rated from the model, and a test section where the major 
portion of the instrumentation is located. Gun bore sizes 
range from 1.1 cm to 3.75 cm (the 2.5 cm gun was used 
in these tests) and the test section is about 1 m across 
and 25 m long. The test section can be evacuated to 
achieve the proper Reynolds number for the model size 
and velocity. In the HFFAF a model of about 50 gm 
mass can be accelerated up to 7 km/s. 


The ballistic range can simultaneously simulate 
the Mach number and Reynolds number of a high-speed 
vehicle in an undisturbed, clean, and well-characterized 
ambient gas. Additionally, there is no sting to effect the 
aerodynamics of the model. As the model flys through 
the test section its position, orientation, and time-of- 
flight are recorded at 16 test stations that consist of or- 
thogonal, focused-shadowgraph systems. In the HFFAF 
the position of the model can be measured to ±0 .006 
cm, its orientation to ±0.1°, and time-of-flight to ±40 
ns. 

In the present study, cones with a bluntness ra- 
tio, dn/db, of .3 and shock generators were launched in 
the ballistic range for the purpose of measuring the aero- 
dynamic coefficients and shock shapes. A drawing of 
the model is shown in Fig. 2. Fig. 3 is a photograph 
of the model and its four piece sabot and pusher plug. 
The sabot protects the model during the high-g’s en- 
countered during launch. After launch, aerodynamic 
forces separate the sabot pieces and the model proceeds 
downrange in free flight The nose of all models was 
fabricated of tantalum to avoid nose-tip ablation and the 
afterbody was fabricated of aluminum. The test condi- 
tions for the experimental data are listed in Table 1. The 
runs were divided into three groups based on freestream 
conditions. In case 1, V = 5 km/ sec. and Rei = 10 5 ; 
in case 2, V = 5 km/ sec and Rei = 10 6 ; in case 3, 

V = 6 km/s and Rei — 10 5 . The test gas was air. For 
all conditions tested, the flow over the models remained 
laminar. When test runs were made at V = 6 km/s and 
Rei = 10 6 , the nose tip temperature became so high 
that the tantalum nose began to spall. These runs are not 
included in this report. 

A new data reduction routine is used to fit the 
calculated trajectory to the measured trajectory using a 
least-squares analysis. The routine assumes a 5 degree- 
of-freedom model of the equations of motion (roll rate is 
assumed constant), written in terms of the aerodynamic 
coefficients. The aerodynamic coefficients are approxi- 
mated by Taylor series expansion in sin <r: 

Cm - (C Ml + Cm 3 sin 2 cr + Cm s sin 4 cr) sin <7 ( 1) 

Cl = (Cli + Cl 3 sin 2 cr + Cl s sin 4 cr) sin a (2) 

Cd — Cd 0 + Cd 2 sin 2 a (3) 

The routine can fit raw data from several test runs si- 
multaneously improving the accuracy of the method. 

The equations and method used are described in more 
detail in the Appendix. 

The data reduction routine provides estimates of 
the error, e, in the aerodynamic coefficients; for exam- 
ple. 




The error estimates take into account the uncertainty in 
measuring model trajectory and freestream conditions 
and the quality of the least-squares fit. However, they 
are sensitive to the number of terms taken in the Taylor 
series expansions used to approximate the aerodynamic 
coefficients. During data reduction, the number of terms 
in the expansions is varied and the set which produces 
the smallest estimated error is reported. The minimum 
estimated error occurs in the region where the experi- 
mental data is concentrated. For the present tests, the 
data is concentrated at resultant angles between 5° and 
10 °. 


FLOWFIELP COMPUTATIONS 

Two numerical algorithms were used to obtain 
the computed results. 13,14 The baseline conditions for 
these computations are listed in Table 2. The first algo- 
rithm, TUFF, is a time-marching scheme and was used 
to compute the axisymmetric, blunt-body results from 
the stagnation point to 70° on the spherical nose. The 
second algorithm, STUFF, is a space-marching scheme 
and was used to compute the remaining afterbody flow 
field. Both codes employ a finite-volume philosophy to 
ensure that the algorithms (including boundary condi- 
tions) are fully conservative. Further, they obtain up- 
wind inviscid fluxes by employing a Riemann solver 
that fully accounts for the gas model used. Either a 
ideal-, equilibrium-, or nonequilibrium-chemistry model 
can be used. 14 Thermal equilibrium is assumed in all 
gas models. The algorithms are total variation dimin- 
ishing (TVD) which allows the computation of flow- 
fields with strong discontinuities without introducing 
any spurious oscillations. They employ a strong cou- 
pling between the fluid dynamic and species conserva- 
tion equations and are made fully implicit to eliminate 
the step size restriction of explicit schemes. Finally, the 
Vigneron approximation is used in the STUFF algorithm 
to allow stable space-marching in the presence of a sub- 
sonic viscous layer. 15 

This combination of a time- and space-marching 
algorithms was chosen since it allowed the solution to 
be obtained on a rather fine grid in a reasonable amount 
of computer time. The blunt body grid measured 29 
cells from the stagnation point to the terminator and 49 
cells of exponentially increasing size in the normal di- 
rection. The first cell was 0 .25 x 10 " 6 m off the body 
surface for the lower Reynolds number computations 
and was decreased by a factor of three for the higher 
Reynolds number cases. The after-body grid used the 
same normal spacing as the blunt body grid. The three- 
dimensional, space-marching computations were per- 
formed with 29 meridonal planes of cells. Nonequilib- 
rium computations at zero angle-of-attack required about 
1.4 CPU-hours on the Cray Y-MP, while computations at 


5° angle-of-attack required over 6.1 CPU-hours. 

The surface temperature as a function of flight- 
time was computed with a one-dimensional, heat-transfer 
model. The surface heat flux was taken to be 


q s (t f ) = 




y/TTtfCXT 


(5) 


An iteration was performed to obtain the wall tempera- 
ture, T 9 , along the model surface and the time-of-flight, 
t /, was taken to be 5 msec. The initial temperature of 
the cone, 7 1 ,-, was assumed to be room temperature and 
the cone’s multiple composition was modeled. A fully 
catalytic wall boundary condition was used. Wall tem- 
perature profiles show that, except at the nose, the wall 
temperature is approximately that of the free-stream. 

RESULTS AND DISCUSSION 

The experimentally determined coefficients of 
the aerodynamic functions are tabulated in Table 3. The 
aerodynamic coefficients at a specific angle-of-attack can 
be obtained from equations 1, 2, and 3. Included with 
the values for the aerodynamic coefficient are values for 
the estimated error. Computations were performed at 0 
angle-of-attack for all cases and at 5 ° angle-of-attack for 
cases 1 and 2. For case 2, three gas models were used. 
The computed aerodynamic coefficients at 0° and 5° 
angle-of-attack are tabulated in Table 4. The experimen- 
tal and computational values are also presented graphi- 
cally for comparison. In these figures, the experimental 
data are plotted as lines and the computed results are 
plotted as points. The number of terms included in the 
expansions of the aerodynamic coefficients and the do- 
main in which these expansions are valid are determined 
by the total angular range of the data. For case 2, the 
angular range extended to a - 8°, and the data is plot- 
ted only over this range. For cases 1 and 3 the angular 
range extended to a - 15°. 

Computed and experimental shock shapes are 
in reasonable agreement. Fig. 4 shows a comparison 
of the computed and experimental outer shock struc- 
ture for case 2. The computed shock is the white line 
superimposed on a shadowgraph from the range test. 

An outline of the body shape used in the computation 
is also included. Unfortunately, the ambient densities 
at the conditions of the tests were too low to produce 
shadowgraphs of sufficient quality to discern the inter- 
nal shock structure produced by the shock generators. 

As will become apparent in the discussion that follows, 
however, the shock generators did have some effects on 
the aerodynamics of the test model. 

Drag 

Experimental and computed drag coefficient 
versus resultant angle-of- attack, cr, is plotted in Fig. 5 



for all three cases and are tabulated in Tables 3 and 4. 
The pressure at the base was assumed to be a vacuum 
for the computational results. Actually, the base pressure 
is estimated to be about half that of ambient, resulting 
in an increase in the computed drag coefficient of less 
than 3%. This correction to the computational result is 
approximately equal to the estimated experimental er- 
ror. Values of Cd calculated at 0° and 5° using mod- 
ified Newtonian theory are also included for compari- 
son. Newtonian theory considers only pressure drag and 
should provide a value close to but slightly below the 
higher Reynolds number case (case 2) as is observed. 
The drag for case 2 is due mostly to the pressure drag 
generated by the blunt nose. The drag due to viscous 
effects is much more important in cases 1 and 3, where 
the Reynolds number is 10 5 . From case 1 to case 2 the 
Reynolds number has increased from Rei = 10 5 to 
Re t - 10 6 while the velocity has remained constant at 
5 km/s. Since skin friction is proportional to Reynolds 
number to the -0.5 power, the viscous drag should be 
three times greater in case 1 than in case 2. This effect 
can be seen in the shift of the experimental drag curve. 
For example, C Do , equal to 0.170 in case 1, decreased 
to 0.112 in case 2(see Table 3). This shift is reflected in 
the computed nonequilibrium-chemistry drag coefficient 
which was 0.1617 for case 1 and 0.1241 for case 2. Lit- 
tle real-gas effects are observed in the computed values 
of Co- This is in agreement with earlier data for blunt 
cones. 

The velocity has increased from V = 5 km/s in 
case 1 to V = 6 km/s in case 2 while the Reynolds num- 
ber has remained constant at 10 5 . There is negligible 
effect of this velocity increase as can be seen in Table 3, 
where the differences in Cd 0 311(1 Co 2 are less than the 
sum of their estimated errors. 

Computational results are in generally good 
agreement with experimental results. Computed drag 
for cases 1 and 3 fall slightly below but within the esti- 
mated error of the experimental values. The computed 
value for Co 0 in case 2 is higher than the experimen- 
tal value by about 10%. This is outside the estimated 
error for that case. The shape of the drag curve with 
respect to angle-of-attack is caused by drag due to lift, 
thus, Co z is expected to be proportional to Cl. Tables 
3 and 4 show that experimental values for Cdz vir- 
tually the same for cases 1 and 2 as are values for Cl 
in the computed results. The computed value of Coz for 
case 1 is 3.46, in good agreement with the experimen- 
tal data; the computed value for Co 2 for case 2 is 2.42, 
which is not in agreement with the experimental data 
The reasons for this discrepancy are presently unknown. 

Moment 

Experimental and computed values of the mo- 


ment coefficient are plotted in Figure 6 with the values 
listed in Tables 3 and 4, respectively. The experimen- 
tal data shows much non-linearity in agreement with 
the results reported by Welsh et al. 7 and Malcolm and 
Rakich 9 for blunt cones with different cone half-angles. 
The non-linearity can also be seen in the values for Cm 
listed in Table 3. For cases 1 and 3, three moment co- 
efficient terms were needed to adequately fit the model 
trajectory. In case 2 only two terms were needed to de- 
scribe the model trajectory due to a smaller range of 
angular motion. Table 3 shows little difference between 
the values for Cm\ 311(1 Cm 2 over the measured range 
of motion for cases 1 and 2 showing little Reynolds 
number effect. This is reflected in the nonequilibrium- 
air values of Cm computed at 5° reported in Table 4. 
While the computed viscous contribution to the moment 
coefficient decreased markedly from 0.0090 in case 1 
to 0.0028 in case 2, the pressure contribution increased 
from 0.0339 in case 1 to 0.0375 in case 2, offsetting the 
decrease in viscous effects. Table 4 shows a difference 
in the moment coefficient in case 2 of 32% between the 
computed ideal-gas value to that of the nonequilibrium- 
air value, demonstrating that the real-gas effects for this 
case made the model more stable. Welsch et al. found 
that while real-gas effects made models with a blunt- 
ness ratio of less than 0.2 unstable, models with a blunt- 
ness ratio greater than 0.2 were stabilized. 7 This agrees 
with the current computational result The calculations 
for case 2 show little change in the viscous contribu- 
tion to moment coefficient between the nonequilibrium- 
air model (-0.0028) and the ideal-gas model (-0.0024). 
The pressure contribution changes from -0.0375 for the 
nonequilibrium air case to -0.0249 for the ideal gas case, 
a change of 66%. The estimated errors for experimental 
moment coefficient for case 3 are too large to attribute 
any effect to the change in velocity. 

Figs. 7 a through c show computed mass frac- 
tions of atomic oxygen for the three cases, respectively, 
at three locations measured from the nose and at 0° 
angle-of-attack. The mass fraction of oxygen in the am- 
bient stream is 26%. The oxygen is fully dissociated 
in the stagnation region of all of the cases, and the fig- 
ures show that significant amounts of atomic oxygen are 
present in the flow along the after body where it has the 
most effect on the aerodynamics. The mass fraction of 
atomic oxygen for all cases decreases to zero at the wall 
reflecting the wall boundary conditions. In case 1 the 
mass fraction is lower than in the other two cases and 
decreases quickly downstream reflecting the lower tem- 
peratures generally present in the flow in that case. The 
calculations for case 2 show little change in the viscous 
contribution to the moment between the nonequilibrium 
air chemistry model (-0.0028) and the ideal-gas chem- 
istry model (-0.0024). The pressure contribution to the 



moment, on the other hand, changes from -0.0375 for 
the nonequilibrium air model to -0.0249 for the ideal 
gas model, a change of 66%. This suggests that real-gas 
effects are important in the pressure field rather than in 
the viscous boundary layer as one might expect. 

Fig. 8 shows simulated oilfiow patterns for 
cases 1 and 2 at a = 5°. The simulated oilfiow pat- 
terns were generated by taking the velocity vector of the 
first gridpoint off the surface. These figures show that 
most of the dissociated fluid from the nose of the body 
is swept to the lee side of the body concentrating atomic 
oxygen on the lee side. Fluid in the boundary layer 
on the windward surface has traversed a weak oblique 
shock and will contain little dissociated oxygen while 
fluid on the lee surface will contain a large amount of 
dissociated oxygen. This partially explains the differ- 
ence between the moment coefficients calculated using 
the ideal gas and nonequilibrium air chemistry models. 
The computations suggest the start of a separation re- 
gion with reverse flow near the second shock generator 
at a = 5°. While the PNS code used to compute the 
flow on this portion of the body can indicate the pres- 
ence of a separated region, it cannot compute flow de- 
tails in that region or the effects on the aerodynamics. 
This separation region is quite small and will have little 
effect on the drag but may effect the lift and moment 
coefficients somewhat. 

The computed solutions are, of course, steady 
state solutions, and the question arises whether the dy- 
namics of the model in ffee-flight in the ballistic range 
will allow the separation to set up at this angle-of-attack. 
The models went through between one and two cycles 
of motion during their 5 msec flight down the range. By 
contrast, the flow time (the time fluid at the nose takes 
to reach the tail) is on the order of 6 p sec. This corre- 
sponds to a reduced frequency of about 0.002. There 
will be enough time for the separation to set up and, 
hence, the trajectory and experimentally determined 
aerodynamic coefficients will reflect the effects of the 
separation. 

Lift 

Figure 9 shows the experimental and computed 
values for Cl for case 2. The model motion is less sen- 
sitive to the lift coefficient than to the other coefficients 
reported in this paper, and in ballistic range testing 1.5 
to 2 wavelengths of motion is required in a given run to 
determine the lift coefficient About half a wavelength 
of motion was obtained for cases 1 and 3. While this 
was enough for the determination of drag and moment, 
it did not allow sufficient confidence in values for the 
lift coefficient. Therefore, lift coefficients are reported 
only for case 2. Computed values of Cl for all cases 
are reported in Table 4. Since our confidence in the ex- 


perimental lift coefficients for cases 1 and 3 was low, 
the effect of variations in the lift coefficient on the other 
aerodynamic coefficients was investigated. A linear lift 
coefficient was assumed, i.e. Cl 3 and Cl s were set to 
zero while Cl { and the drag and moment coefficients 
were allowed to vary. Values of Cl equal to 1.5 times 
the computed value were input into the reduction routine 
and the moment and drag coefficients were recalculated. 
This resulted in changes in Cd and Cm on the order 
of 0.05%. These variations are much smaller than the 
expected accuracy of the drag and moment coefficients 
themselves. 

Agreement between the computed and experi- 
mental values of lift coefficient at 5° for case 2 are ac- 
ceptable. A value obtained using modified Newtonian 
theory is included for comparison. This value is very 
high as is expected. The computed nonequilibrium air 
value of 0.0927 is greater than the ideal gas value of 
0.0778 showing a significant real-gas effect. The real- 
gas effects for the 5° cones are much greater than that 
reported by Welsh et al. for 10° cones. 

The calculated value of Cl for case 1 (see Ta- 
ble 4) was 0.0908, only 2% below the value for case 2. 
This suggests that Reynolds number effects are not as 
important in the lift coefficient as was suspected. 

Damping 

Damping is due to 1) a slightly different angle- 
of-attack from forebody to afterbody due to the curved 
flight path, and 2) a time lag before a translation in nose 
position can result in a translation of the shear flow 
on the afterbody. The first order damping coefficient, 
CM d = Cm, + CM,* is extremely difficult to obtain in 
most facilities. Since all of the aerodynamics are re- 
flected in the motion of the model, the damping coef- 
ficient must be included in the data reduction routine. 
When there are onlu a few cycles of motion, the trajec- 
tory of the model is not very sensitive to the damping 
coefficient, thus, the estimated errors obtained for these 
coefficient are typically large. However, designers of- 
ten need only know the sign and order of magnitude 
of these coefficients to properly design a vehicle. The 
damping coefficient obtained for all three cases are in- 
cluded in Table 3. Note the large differences in damping 
coefficient between the three cases. Unsteady codes are 
needed to compute the damping coefficient; hence, none 
were computed for this model. 

CONCLUSIONS 

Experiments and computations have been per- 
formed on 5° half-angle cones with a bluntness ratio 
of 0.3 and with shock generators at 5 and 6 km/s and 
at Reynolds numbers of 10 5 and 10 6 in air. Experi- 
mentally derived aerodynamic coefficients have been 
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determined for drag, moment, lift, and damping with 
estimated errors. The computations employed ideal-, 
equilibrium-, and nonequilibrium -gas chemistry mod- 
els for air and realistic boundary conditions at the wall. 
The calculated moment and lift coefficients demonstrate 
real-gas effects. These effects increase the stability of 
this test model, in agreement with previous results. The 
strong shock produced by the at these speeds generates 
temperatures high enough to dissociate all of the oxy- 
gen in the blunt nose region. The computations show 
that significant amounts of atomic oxygen persist in the 
flow along the after body where it can effect the aero- 
dynamic performance of the model. The presence of 
the atomic oxygen is seen to have a large effect on the 
pressure contribution to the moment coefficient and little 
effect on the viscous contribution. The experimental re- 
sults have been compared with the computations. With 
few exceptions, there is good agreement between exper- 
imental and computed shock shapes and aerodynamic 
coefficients. The computations can provide detailed in- 
formation about the state of the gas. While measurement 
techniques have been developed which can provide ex- 
perimental confirmation of these details, they have only 
recently begun to be applied in hypervelocity facilities. 

APPENDIX 

A Five-Degree of Freedom Routine 
for Determining Aerodynamic Coefficients 
from Free-Flight Experiments 

Aerodynamic coefficients can be obtained from 
free-flight experiments with a least-squares procedure 
that fits functions describing the motion of a projectile 
to experimental data. 16 The procedure has been modified 
to a weighted least-squares procedure. The motion of 
the model is described by twelve first-order differential 
equations. 17 For a five-degree of freedom system, the 
roll rate is assumed constant and the system reduces to 
eleven first-order differential equations: 
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+ coriolis and gravity terms ( 17) 

L 

q w = — — + conolis and gravity terms (18) 

mV 

Y 

r w = — + conolis and gravity terms ( 19) 

mV 

and ( P \jj , Qw i Ru>') — ( Pw > Qw > ri^/) plus additional terms 
due to a rotating earth. In these equations, V is the ve- 
locity of the projectile, (x, y,z) is the location of the 
projectile in reference to an earth fixed coordinate sys- 
tem, of and are the pitch and yaw angles relative to 
the wind axes, p, q, and r are the roll, pitch, and yaw 
rates, and & wt and V'w give the orientation of the 
wind axes to the earth fixed axes. The model specific 
parameters m, I y , and I x are, respectively, the mass, the 
moment of inertia about a transverse axis through the 
center of gravity, and the moment of inertia about the 
axis of symmetry. The aerodynamic information is con- 
tained in the functions describing the pitching and yaw- 
ing moments, M and N, the functions describing the lift 
and yawing forces, L and Y , and the function describing 
the drag, D. 

For an axially symmetric projectile, these aero- 
dynamic functions can be approximated by the expan- 
sions: 


M = 2 (Cm 0 + Cm, cos (3 sin a) 


+ CM,q + Cm, a 

(20) 

N = + CU ' sin 0) 


+ Cm,t COS a - Cm*P 

(21) 

L- 2 (C Lo + C L , cos (3 sin a) 

(22) 

K=£^(C ro + C i ,sin/3) 

(23) 
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D = 2 C ° 

(24) 


where p is the density of the fluid, A is the base area 
of the model, and d is the diameter of the base. The 
coefficients Cm 0 > Ch 0 * Cl 0 > and Cy 0 are the moments 
and lifts at zero angle-of-attack and are caused by small 
asymmetries in machining. The aerodynamic coefficients 
Cm , > Cl,, and Cd are assumed to be functions of a 
(the resultant angle-of-attack), the Mach number, and 
the Reynolds number. The aerodynamic damping terms, 
Cm, and Cm, are assumed to be constant. For the test 
conditions reported in this paper, j3 % — rcos a and 
6> m q. It is therefore impossible to determine Cm, and 



Cm & separately; only the sum of these two coefficients 
can be determined. The other aerodynamic coefficients 
are approximated by 

Cm* = Cmi + Cm, sin 2 <7 + Cm s sin 4 cr (25) 
Cl* = C L , + Cl, sin 2 a + Cl , sin 4 <7 (26) 

Co = Coo + Cpi s i° 2 & (27) 

The estimated errors for the aerodynamic coef- 
ficients and for functions of these aerodynamic coeffi- 
cients are readily available from the least squares analy- 
sis. For instance, if the angle-of-attack is known exactly, 
the variance for Cd is equal to 

Var ( C D ) = Var(Co 0 ) + 2Cov {Cd^C^) sin 2 ct 


+ Var ( Cd 2 ) sin 4 a (28) 

where the variances and covariance are given by 

Var(C Ds ) = ^ci„.Cr„ Var ( I ) < 29 > 

Var (Co,) = , Cdj Var(x) (30) 
Cov (Co 0 ,Co,) = ^5' j , CDj Var(i) (31) 


The variance and, hence, the predicted error for the drag 
are functions of the variances for each coefficient, their 
covariance, and the angle-of-attack. Since the covari- 
ances are often negative, the minimum predicted error 
can occur at non-zero angles-of-attack; it typically oc- 
curs in the region where the angular data is clustered. 
Similar results can be obtained for both the lift and mo- 
ment coefficients. 
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Tabic I. Experimental Conditions 



Shot # 

V (km/s) 

Rej 

p(kp B ) 

Tempt K ) 

Mass(gm) 

Length(cm) 

Base diam(cm) 

■wmi 

Case l 


4.99 


0.878 

297.4 

Hi 

3.62 

0.886 

sis 


B 

4.98 

B nriCTi 

0.884 

SB 
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3.62 

0.886 
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9.443 
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5.00 

IKYCTTiM 
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3.63 


1.85 

Case 3 


6.08 

iffPIFTTOBF 

0.793 

294.2 

2.09 

3.64 

0.889 

1.88 

MppaEp; 


6.05 

lytfirfTMB 

0.800 

294.2 

2.08 

3.64 

0.889 

1.88 


Tabic 2. Baseline Conditions for Computations 



V (km/s) 


Temp(K) 


Case 1 

5.0 

0.0103 

298 

1.813 

Case 2 

5.0 

0.1126 

298 

1.813 

Case 3 

6.0 

0.0103 

298 

1.813 


Table 3. Experimental Aerodynamic Coefficients 



V (km/s) 


Rej 

Cd 0 


Cl, 

Clj 

Cmi 

c Mi 

C Mi 

Cm/ 

Case 1 

5 

14.4 

10 5 

0.1700 

±0.0065 

3.532 

±0.184 



-0.204 

±0.028 

-30.7 

±1.7 

309. 

±20. 

-2.24 

±0.36 

Case 2 

5 

14.5 

10 6 

0.1124 
±0.001 1 

3.537 

±0.140 


38.0 

£4 .6 

-0.243 

±0.017 

-30.0 

±1-6 


-0.30 

±0.04 

Case 3 

6 

17.6 



10 s 

0.1627 

±0.0104 

3.269 

±0.152 





410. 

±54. 

-11.45 

±2.00 


* C m 4 — Cm, + C Mi 


Table 4. Computed Aerodynamic Coefficients 




Rej 

cr(deg) 

gas model 


C L 

Cm 

Case 1 

5 

mm 

0 

non-eq. air 

0.1617 
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5 

non-cq. air 

0.1880 

0.0908 
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0 

non-eq. air 

0.1241 



■ 



5 

non-eq. air 

0.1425 

0.0927 

-0.0404 





ideal gas 

0.1272 






5 

ideal gas 

0.1413 

0.0778 

-0.0275 




0 

eq. air 

0.1263 



Case 3 

6 

HUSH 

0 

non-eq. air 

0.1592 









































































































































DRAG COEFFICIENT 
CASE 



1, EXPERIMENT 

□ 1, NON-EQUILIBRIUM GAS 

2, EXPERIMENT 

O 2, NON-EQUILIBRIUM GAS 
A 2, PERFECT GAS 
V 2, EQUILIBRIUM GAS 



o. ANGLE OF-ATTACK, deg 


3. Photograph of model with sabot pieces. - 5. Experimental and computed Co vs angle-of-attack. 




4. Shadowgraph of 5° model in free-flight at V — 5km/s 6. Experimental and computed Cm vs angle-of-attack. 

and Re = 10 4 * 6 (case 2). The white lines superimposed 

on the photograph are the computed shock shape and 
body outline. 




NORMAL DISTANCE, 




0 2.5 5.0 7.5 10.0 

a, ANGLE-OF-ATTACK, deg 


9. Experimental and computed Cl vs angle-of-attack for 
case 2 (V = 5 km/ s; Re = 10 6 ). 
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for the 3-D Navier- Stokes Equations 
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Abstract 

One method of solving the Navier* Stokes equations 
about complex and realistic aerodynamic configurations is 
to use a zonal method. In this method the overall flow 
field domain is subdivided into smaller blocks or zones. In 
each of these zones, the flow field is solved separately of 
the other zones. The boundary data for each zone is pro- 
vided by the neighboring zones. The major difficulty of 
the zonal methods has been how to maintain overall con- 
servation for arbitrarily shaped zones. A new method of 
conservative patched zones has been developed. It uses 
structured meshes in the individual zones. The interface 
between the zonal block faces is defined by the union of 
the face points of adjoining blocks. An unstructured grid 
is generated upon which the interface fluxes can be de- 
termined. Flux balancing of the interface fluxes is then 
easily achieved to obtain global conservation. The method 
has been implemented into two Navier-Stokes codes. The 
use of the procedure is easily implemented into other fi- 
nite volume codes. There are no topological restrictions 
on the zonal boundaries; e.g. the zonal interfaces can be 
curved surfaces for ease m constructing structured meshes 
in each of the zones. Several examples are presented to 
demonstrate the viability of the interfacing procedure. 

Introduction 

There are two basic approaches of numerically simu- 
lating the Navier-Stokes equations about complex and re- 
alistic aerodynamic configurations. One is based on struc- 
tured meshes in which the neighbors of a mesh point are 
known implicitly. The other approach is based on an un- 
structured mesh in which the neighbors of a mesh point are 
not known implicitly and this information must be stored 
for each point. Numerical methods based on structured 
meshes are well developed, but suffer limitations when 
dealing with complex configurations in that it is difficult, 
if not impossible, to generate a single mesh about such a 
configuration and still have the required mesh qualities for 
stable and accurate numerical solutions. The generation 
of unstructured grids about complex configurations is, in 
principle, much easier; however, the numerical methods for 
such grids are not yet mature enough to compete against 
structured grid numerical methods. 
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The structured grid generation problem for complex 
configurations can be alleviated by zonal methods. In 
this method the overall domain is subdivided into smaller 
blocks or zones. In each of these zones, a grid is gener- 
ated and the flow field is solved independently of the other 
zones. The boundary data for each zone is provided by 
the neighboring zones. The major difficulty of the zonal 
methods has been how to maintain overall conservation for 
arbitrarily shaped zones. The resolution of this difficulty 
is the purpose of this paper. 

There are two types of zonal methods in common us- 
age today. One is the overlaid zones in which both zones 
share a Common interface region, as used in the Chimera 
[1] approach. The other is the patched zonal technique 
where the two zones share only a common boundary. Ex- 
amples of this approach are the zonal method of Rai [2] and 
Thomas [3].- Both methods have advantages and disadvan- 
tages depending on the application. With patched grids it 
is much easier to maintain conservation since there is only 
one boundary across which the two zones communicate. 
On the other hand for moving multiple bodies the overlaid 
grid approach has certain advantages if conservation is not 
important, i.e., the flow field is continuous with no shock 
waves or shear surfaces. If conservation is important, as it 
is in this investigation, then we are restricted to patched 
zones. 

Most methods using conservative patching methods 
are restricted to interface surfaces which are planar sur- 
faces due to the minor gaps and overlaps that occur at a 
curved interface if the two zones are not mesh continuous. 
Furukawa et al [4] attempted to resolve this problem by 
using only one of the zones to determine the zonal bound- 
ary for both zones. The open question that remains is 
which zone determines the interface boundary. Furukawa 
et al chose to use the zone which has the better resolution. 
This can result in loss of accuracy if the mesh ratio changes 
in the interface, for example at a viscous boundary layer. 

In this paper the zonal interface boundary is deter- 
mined by the union of all the face mesh points of both ad- 
joining zones. The interface surface is now unique and de- 
termined much more accurately than either one of the indi- 
vidual face surfaces. The collection of interface points is in 
general no longer structured and readily available unstruc- 
tured grid generation techniques can be used to construct 
(triangulate) an interface grid. With the triangulated in- 
terface grid, the metrics (i.e., surface area normals) and 
cell volumes can be determined for each of the interface 
cells. 

The development of the interface algorithm is dis- 
cussed in the following sections. The procedure haS 
been implemented into two finite volume three-dimensional 
Navier-Stokes codes, namely TUFF [5] and a finite volume 
version of CNS [6] henceforth called CNSFV. We will limit 
the exposition to the diagonal version of the Beam and 
Warming scheme [7] as used in the CNSFV code, but re- 
sults from the TUFF code will also be presented. 
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Numerical Scheme 
Navier- Stokes Equations 


The three-dimensional thin-layer Navier-Stokes equa- 
tions in strong conservation law form in curvilinear coor- 
dinates are 

drVQ + dtF + drjG + dcH = Re-^S (1) 

where 
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The contravariant velocity components are defined as 

U = V(£ t + CxU -f £ y v + £ z u>), 

V = V(r) t + T) x u + 7} y v + r] 2 w), 
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and the viscous flux is given by 
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( 4 ) 


The metrics used above have a different meaning for a 
finite volume formulation compared to the finite difference 
formulation of [ 6 ]. Referring to a typical finite volume cell 
as shown in figure 1 , the finite volume metrics are defined 
(see, for example, Vinokur [ 8 ]) as 


s ;+i = + 


-[(r 7 - r 4 ) x (r 8 - r 4 ) 4 - (r 3 - r 4 ) x (r 7 - r 4 )] 


s fc+| = -Si.fc+i* + *y,k+±j + 


-[(r 7 - r 2 ) x (r 3 - r 2 ) +'(t 6 -r 2 ) x (r 7 - r 2 )] 


s z+i — + 5 »,j+ + • s *,i+i k ( 7 ) 


= 2 [( r " - r s) x (r 6 - r 5 ) 4- (r 8 - r 5 ) x (r- - r 5 )] 


The finite volume metrics represent the cell face area 
normals in each of the curvilinear coordinates (£, 77 , £). 
They are related to the metrics introduced in equations 
(1 - 5) as follows 


CxV = s x,j+% 
iyV = S y j + 1 


with 


m l = £i 2 + £y 2 + £z 2 , 

m 2 = C x u ( 4- CyV< 4- ( z w<, (5) 

m 3 = ( u 2 4 - v 2 4 - w 2 )/ 2 4- (Pr( 7 - l)) - 1 (a 2 )^ 


The pressure is given by the equation of state 

p = (7 - l)(e ~ pi? 2 + v 2 4- w 2 )/2) 
Metric Terms 


( 6 ) 


VzV = S t , k+ 1 
T]yV = <9 y ,Jk + i 
VzV = s Ztk+ 1 

C.^ = v*+J 
CyV = -Sy^+I 
<*V = 


( 8 ) 


2 



The volume of the computational cell is given by 


6qG = B = TjjAfjTj, 


V = l [(r 4 - ri) x (r 3 - ri) ■ (r 7 - r 8 ) 
o 

+(r 8 “ ri) x (r 4 - ri) • (r 7 - r 8 ) 

+(r s - r a ) x (r 2 - r x ) * (r 7 - r 8 ) 

+(r 2 - ri) x (r 6 - r x ) • (r 7 - r e ) 

+(r 6 - ri) x (r 8 - ri) • (r 7 - r 6 ) 

+(r B - r x ) x (r 8 - ri) - (r 7 - r 8 )] (9) 

and is the finite volume equivalent of the inverse Jacobian 
of the coordinate transformation in the finite difference 
formulation of [6]. 

Diagonal Beam- Warming Algorithm 

The implicit Beam- Warming scheme for the finite 
volume formulation is given by 


d Q H = C = IVA^- 1 

The T(,T n ,T^ are the eigenvector matrices of A, B, C, 
respectively with as the respective eigenvalues. 

We also have 


N = T{~ l T v 
P = T n ~ l T( 

Each of the factors of the implicit operator of equa- 
tion (12) has an artificial term added to stabilize the cen- 
tral difference operator. The dissipation is based on Jame- 
son’s nonlinear second and fourth order dissipation and for 
the operator takes the form 


hV - 1 Vi [*j+ 1(* (2) A* • - e (4) A* Vi A* -)]AQ n (13) 
with 


e (2) = K 2 max(-fj+i (I 4 ) 


[/- V-'h(6tA n + 8 v B n + 6 (: C n -Re- 1 6 < M n )}6Q-~ i = R n 

( 10 ) 


where 


_ \pj + 1 - + Pi-i i 

|pj+i + 2 Pi + Pi- 1 1 


(15) 


R n = -V- l h[6 i F n + S v G n + S^E 71 - Re~H ( S r -} (11) 

The convective three-dimensional flux Jacobians A, 
B. C and the viscous flux Jacobian are defined in the ap- 
pendix of [5]. With the use of approximate factorization 
and diagonalization of the flux Jacobian matrices, a scalar 
pentadiagonal algorithm [7] can be derived as 


r € [J + v~ l hit A«]JV[/ + V-'W, A,]P. 

[i + y- 1 w < A c ]r c “ 1 AQ n = i? n (12) 

where S ^ is a central difference operator and A Q n = 

Q n+1 — Q n with Q n+1 — Q(f n + A). The viscous terms 
are not included in the left-hand (implicit) side. The ar- 
tificial dissipation is included in both sides and is derived 
below. 

The inviscid flux Jacobians are diagonalized as fol- 
lows: 


= max(0,/C4 — e^) (16) 

where K 4 are constants of o(l), and V{ 3X6 the f° r_ 
ward and backward difference operators. is a modi- 

fied spectral radius defined els 

d j+ i = Oj + Vj +i 

Cj = <Tj( 1 + y] max{a k / a j , <x,/<Tj )) , (17) 

<tj = [\U\+ay/s x i + s v *+s z *]' (18) 

and where cell centered surface areas are used, e.g. 


d Q F = A = T i A i T i ~ l 


3 


The modified form of the spectral radius, equation 
(17), is suggested by Turkel [9] to account for large aspect 



ratio computational cells as for example in a viscous layer. 
Similar dissipation terms are obtained for the rj— and £ — 
operators. The dissipation terms added to the right hand 
side of equation (10) are identical to those given above 
except that A Q n is replaced by Q n . 

To show that the above scheme is conservative and 
has the telescoping property, the right hand side of (12) 
can be written as 

\ 

R n = -hV-'lF”,! 


+ff" , - ff," • - Re~'S? . + R‘~'S? ,1 (19) 

where F.G.H.S are the numerical fluxes and are defined, 
for example for F, as t . 


J = + 5 ;+ i(‘ <!) A (Q n - e <4) A £ Vs A { e») i+1 

■V, (20) 

' V 

with 


F? +i = j [ s i,j+ + Si) " + a »j+}(Sj+i +SiT 


( 21 ) 

Similar terms are obtained for the numerical fluxes 
in the other two coordinate directions. 


Boundary Conditions 

To complete the equation set. boundary conditions 
must be specified. With the use of curvilinear coordinates, 
the physical boundaries have been mapped into compu- 
tational boundaries, which simplifies the application of 
boundary conditions. The boundary conditions to be im- 
plemented for external viscous or inviscid flows include (1) 
inflow or far field, (2) outflow, (3) inviscid and (4) viscous 
impermeable wall, and (5) symmetry conditions. For ex- 
ternal three-dimensional flow fields about closed bodies, 
the topology of the grid usually introduces (6) grid singu- 
larities which require special boundary conditions. The use 
of zonal methods can avoid the generation of grid singu- 
larities, but requires (7) special zonal interface boundary 
conditions. For compressible flows these zonal boundary 
conditions should be conservative to maintain global con- 
servation. 

In the finite volume approach, the specification of 
boundary conditions reduces to specifying the appropriate 
numerical fluxes at the boundaries. The details of im- 
plementing boundary conditions (1) through (5) are well 
known and are given in [5] and [10]. The grid singularity 
boundary condition is described below and the interface 
boundary conditions are given in the next section. 


The grid singularity boundary condition is similar to 
the symmetry boundary condition for the inviscid and vis- 
cous fluxes in that there is no flux through that boundary. 
If, however, that is all that is done the results shown in 
figure 2a are obtained. These results are the density con- 
tours of a Mach 8 viscous blunt body flow. As shown in the 
figure, a nonphysical behavior appears at the singular line. 
The nonphysical results are due to a local violation of the 
entropy condition [11]. For central difference schemes, the 
artificial dissipation is the only stabilizing (entropy pro- 
ducing) mechanism available. At the grid singularity, the 
spectral radius, and hence the artificial dissipation, van- 
ishes due to the vanishing of the metrics. The introduction 
of Harten’s entropy correction [12] resolves the difficulty 
and the results are shown in figure 2b. 

Interface Method 

Most methods using conservative patching methods 
are restricted to interface surfaces which are planar due 
to the minor gaps and overlaps that occur at a curved in- 
terface if the two zones are not mesh continuous. This is 
demonstrated in figure 3a. Furukawa et al [4] attempted to 
resolve this problem by using only one of the zones to de- 
termine the zonal boundary for both zones. This is shown 
for the two-dimensional example in figure 3b. The open 
question that remains is which zone determines the inter- 
face boundary. Furukawa et al chose to use the zone with 
the better resolution. This, however, can result in loss of 
accuracy if the mesh ratio changes in the interface as for 
example at a viscous boundary layer. 

In this paper the zonal interface boundary is deter- 
mined by the union of all the face mesh points of both ad- 
joining zones as shown in figure 3c for the two- dimensioned 
case. The interface surface is now unique and determined 
much more accurately than by either one of the individ- 
ual face surfaces. This presumes that the individual points 
from both zones lie in the interface surface and that the 
interface surface is itself smooth and continuous. The col- 
lection of interface points is, in general, no longer struc- 
tured, and readily available unstructured grid generation 
techniques can be used to construct (triangulate) an inter- 
face grid. With the triangulated interface grid, the metrics 
(i.e., surface area normals) and cell volumes can be deter- 
mined for each of the interface cells, i.e., the computational 
cells that touch the interface surface. 

The redefinition of the interface surface has modified 
all the interface cells and they are no longer hexahedral, 
but rather multifacetted. Determining the surface area 
normals and cell vol um es now becomes more complicated. 
There are three types of interface cells, namely, face cells, 
edge cells, and comer cells. The face cells have only one 
face bordering other zones, whereas the edge and corner 
cells have two or three faces in contact with other zones, 
respectively. A typical cell at the interface is shown in 
figure 4. 

The modification introduced by the interface surface 
requires that the metrics and cell volumes of the interface 
cells be corrected to account for the changed shape of the 
interface cells. Corrections are required for the area nor- 
mals of the cell face touching the interface surface, the four 
sidewalls, which may no longer be quadrilaterals, and the 
cell volume. For inviscid steady flows, the cell volume cor- 
rections are not needed since the volume has no effect on 
the steady solution as shown by equation (11). However, 
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for all viscous flows and inviscid unsteady flows the volume 
corrections are necessary. 

There are two ways that the interface grids can be 
triangulated. The first involves eliminating all the grid 
lines from the structured grid cell faces and constructing 
the triangulated mesh. Since in general the new grid lines 
will not be aligned with the original grid lines, a clipping 
algorithm is used to clip those triangles which lie outside a 
particular individual face cell of either zone. An example pf 
this triangulation procedure is shown in figure 5. Here two 
zones have square faces each consisting of a Cartesiaii.^ftd 
of 15 x 15 uniformly spaced points. The two faces are ori- 
ented at an angle of 45° to each other (see figure 5a). The 
union of both sets of face points results in an unstructured 
collection of points which are then triangulated with an 
advancing front unstructured mesh generation procedure 
[13]. The resulting mesh is shown in figure 5b. 

The second procedure retains the original grief lines 
and triangulates any set of points which form a pplygo^ of 
more than three sides. This approach avoids the use. pf a 
clipping algorithm since each face cell contains an ip^egral 
number of interface triangular cells. An example of this 
procedure is shown in figure 10b for two zones with polar 
grids. Valters. 

The unstructured interface grid requires <^ k£t 
of pointers be defined. These pointers indicate whfdh two 
interface cells share a common interface area or section. 
These sections need not be triangular even though the sur- 
face has been triangulated. A common section (or poly- 
gon) may be composed of several triangles if the surface 
is triangulated by the second procedure described above 
or it may be composed of several clipped triangles if the 
first procedure is used. For a more efficient intens.ee flux 
computation it is convenient to define another set c: cross- 
reference pointers. These cross-reference pointers identify 
which interface polygons are in contact with eacn of the 
interface cells. 


Numerical Interface Flux 

The final step required for the conservative interface 
algorithm is to determine the numerical fluxes at the inter- 
face. Originally it was anticipated that the interface flux 
could be determined by any stable numerical scheme, not 
necessarily the one used in the interior of the zone. How- 
ever, it was found that if the numerical schemes differ, then 
” glitches” always appear at the zonal interface. It was also 
necessary to maintain the same order of accuracy for the 
interface scheme as for the interior scheme. An example 
of this case is shown in figure 6. Here the interface flux is 
determined by a first order upwind scheme and the interior 
by a second order upwind scheme. The unfortunate results 
are the discrepancies at the interface boundary as shown 
by the pressure contours. The use of the same second or- 
der upwind scheme for the interface flux eliminates these 
discrepancies. 

The above example shows that it is not possible to 
completely generalize a conservative interface algorithm. 
The numerical flux must be determined by the same nu- 
merical scheme as used in the interior of the zone. The 
geometric aspects of the interface algorithm can be gen- 
eralized, which is essentially the most difficult part of the 
interface procedure. The numerical fluxes at the inter- 
face must, (and should) be computed by the particular flow 


solver involved. 

The computation of the interface numerical flux is 
relatively straightforward. Since it is scheme dependent it 
is given is general form first, valid for many schemes up 
to second-order accuracy. A specific form will be given for 
the CNSFV code. Figure 7 shows two zones, jl and j2, and 
the interface with the individual surface polygons labelled 
by ”i”. The flux at the interface section ”i” is given by 


= $[«* ,gji=i,9ji=2>9j2=i»9j2=2] 

where is the function defining the particular scheme un- 
der consideration, s' is the area normal of the interface 
section ”i”, and the qjj and qj 2 are the cell-centered values 
in zones 1 and 2, respectively. 

The numerical flux at the opposite face of the inter- 
face cell is slightly more complicated. It is 


Fj2= | = $[fj2=f i£»2=l>9j2 = 2, £-2=3, 9jl = l] 

where q jJ= i is the area (i.e., |^|) weighted average of all 
the interface cells of the zone ”1” sections that comprise 
the face cell of zone ”2” touching the interface surface. 
If higher order schemes are considered, then special care 
must also be taken for the fluxes at the next level of cell 
faces, ie. s. The above formulation is valid for all 

cell-centered finite volume schemes. 

For the particular case of the CNSFV code, the in- 
terface numerical fluxes axe 


F? = F? + o’ i (e (2) A* Q n - e (4) A* Vi A { Q n ) i 

3 3 s 3 

where 


F? = ^[s‘x(/j2=i + /ji=i) n + s 'v(9j 2=1 +9ji=i) n 
2 Z 

+S X z {hj2-\ 

and 


O I = Oj\—\ + <Xj2=l 

The flux at j2 = | is similarily determined, taking care to 
use the appropriate area averaged values of all cell-centered 
values. 

Results 

To validate the interface algorithm, several tests were 
conducted. The first test is the freestream preserving test. 
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In this test, the inflow and and permeable wall boundary- 
conditions (ie. conditions (1),{3) or (4) from the boundary 
conditions section) are set to the freestream conditions and 
the solution is converged. If the initial conditions were also 
set to the freestream condition, then the residual should 
be at machine zero (R= O(10 -14 )) and remain there for 
all subsequent iterations provided that the flow field is dis- 
cretized properly with no gaps or overlaps in any of the 
computational cells and interface boundaries. Indeed, for 
the CNSFV code, this was the case. 

However, the freestream test is not a good indication 
of the accuracy of the scheme. Because of the telescoping 
property of the scheme, the surface area normals can be 
computed inaccurately (even erroneously) and the scheme 
will still pass the freestream test. To test for accuracy, 
a single cell residual, equation (19), is computed with the 
same freestream conditions imposed as above. In this case, 
the maximum residual was of the C?(10~ 8 ) on a 64 bit 
machine (Cray YMP). If the grid and the metric terms 
are computed in double precision (ie. 128 bits), then the 
maximum residual reduces to O(10 — 12 ). This indicates 
that machine roundoff errors are not yet a problem, but 
can be if the meshes are refined much further. 

Three different flow cases covering the entire Mach 
regime from incompressible, to supersonic, and to hyper- 
sonic flows with finite rate chemistry were computed with 
the TUFF code. The same basic conservative interface al- 
gorithm described above was used in all three cases, how- 
ever the conservation law equations differed for each of the 
three cases. 

The first case involved an incompressible inviscid flow 
about a cylinder. The two zone mesh is shown in figure 8a. 
The pressure contours are given in figure 8b and the surface 
pressures in 8c. The results across the zonal boundaries are 
smooth and continuous. 

The second set of results are for supersonic blunt 
body flow. This case is an inviscid Mach 2, axi symmetric 
blunt body flow computed on the four-zone mesh shown in 
figure 9a. Figure 9b shows the solution and the bow shock 
position. For these results, the analytic shock location at 
steady state [14] is shown by the solid squares. The solu- 
tion is shown in terms of mach contours on a background 
grid which is cell centered for plotting purposes only. As 
can be seen the multizonal computed and analytic shock 
shapes compare quite well. 

The final flow results obtained are for a viscous hy- 
personic flow about a hemisphere at Moo = 15-3 and 
Re = 2.2 x 10 5 /m. The interface triangulation for the 
two zones containing the grid polar singularity is depicted 
in figure 10b. The flow results in terms of Mach and 
atomic oxygen concentration contours are shown in figures 
10c and d, respectively. Again, the results indicate that 
the solution contours are smooth and continuous across 
the zonal boundaries. Although not shown, the computed 
shock stand-off distance agreed well with the experimental 
data of reference 15. 

Closing Remarks 

A conservative zonal interface algorithm has been 
presented. It uses some of the best features of both the 
structured and unstructured mesh CFD technology. The 
interface surface grid is unstructured from which the met- 
rics and interface fluxes can be readily constructed to ob- 


tain the proper conservative interface algorithm. For effi- 
ciency and rapid convergence, the flow solver within each 
of the zones is based on structured mesh CFD technology. 
The ordering between the zones can be unstructured for 
maximum flexibility in constructing zones and grids about 
complex and arbitrary configurations. The interface algo- 
rithm has been implemented into two three-dimensional 
Navier-Stokes finite volume codes (TUFF and CNSFV) 
and has shown to yield good results. Further testing is 
being conducted for more complex and realistic aerody- 
namic configurations. The procedure is general and can 
be easily implemented into other finite volume codes. 
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a) grid singularity problem. 



b) grid singularity problem resolved 

by Harten’s entropy correction [12]. 


Figure 2. 


Density contours for a viscous blunt body flow 
Moo = 8.0, Re D = 2.19 x 10 5 ,a = 5.0° 




a) "classical” 

works only for planar interfaces, 
otherwise conservation errors occur. 
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b) Furukawa et al [4] method 
uses boundary of finer mesh 
zone to determine the interface boundary. 



Figure 4. A typical zonal interface finite volume cell. 
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Figure 5a. The two structured face grids 
of a typical interface surface. 
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c) present procedure 

uses all the interface points 
to define the interface boundary. 


Figure 3. Zonal Boundary on a Curved Interface 



Figure 5b. The unstructured grid 

of a typical interface surface. 
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Pressure Distribution on a Cylinder 



a) Four zone grid. 


b) Mach contours. 

■ - analytic bow shock location 
(Lyubimov, NASA TT-F715, 1973). 


Figure 9. 


Axisyrnmetric inviscid blunt body flow at Moo — 2.0. 
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(c) Mach Contours 


(d) Atomic Oxygen Concentration Contours 


Figure 10. Viscous Hypersonic Blunt Body Flow with Finite Rate Chemistry M = 15.3, Re = 2.2 x 10* /m. 
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Abstract 

A computational approach for high-speed base 
flows of reentry shapes is demonstrated. The 
approach is based on solving the thin-layer Navier- 
Stokes equations for perfect gas, equilibrium air, or 
chemical nonequilibrium air. An efficient space- 
marching algorithm is used whenever possible on the 
body and in the far wake. A blocked time-marching 
algorithm is used to calculate the embedded 
subsonic flow in the base region. A two-equation 
turbulence mode! which includes compressibility 
effects is used to extend the approach to Reynolds 
numbers above the transition to turbulent flow. 
Comparisons with data are made for both laminar and 
turbulent base flows, with good agreement in all 
cases. 

Introduction 

Base flows and their effect on the wakes of 
reentry vehicles have long been a subject of both 
experimental and theoretical research. The signature 
of a reentry vehicle is to a large 
extent determined by the character 
of its wake, which extends hun- 
dreds of base diameters behind 
the body. The properties of the far 
wake of a slender body at high 
speed are, in turn, strongly depen- 
dent on the near-wake flow field. 

Only limited theoretical models 
currently exist for the flow in the 
near-wake region. Computational 
tools, are needed to aid in the 
analysis of these types of flows. 


The flow field about a slender reentry vehicle can 
be divided into four regions as is illustrated 
schematically in Fig. 1. Because of the strong 
dependency of the near-wake flow field on the flow 
over the body, calculation of the base flow must 
begin with a computation over the body itself. Many 
techniques are available to calculate the flow field 
over the body. Space-marching parabolized Navier- 
Stokes (PNS) codes can be used to efficiently solve 
for the flow field over the body. The intermediate and 
far wake, the most downstream region of Fig. 1, can 
also be efficiently computed using a PNS code, 
provided initial conditions can be specified. The near 
wake, or base region, encompasses the separated 
flow region, the recompression, and the acceleration 
to supersonic flow. This region is difficult to model. A 
time-marching Navier-Stokes (TNS) computation is 
usually required to model the flow field in this area. 
This paper is primarily directed toward base-flow and 
near-wake computations. The overall objective of the 
effort, however, is directed toward prediction of flow 
quantities, chemical species, and electron concen- 
trations in the intermediate and far wake. 
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High-speed base flows are extremely rich in fluid 
dynamic and thermochemical complexities. Separa- 
tion and reverse flow are the principal features in the 
base region that complicate flow-field predictions. 
The local Mach number may vary from above 10 to 
less than one across the shear layer. Strong shocks 
and expansions produce large temperature gradients 
in the flow. Temperatures which are high enough to 
produce dissociation and ionization in air occur. If 
transition to turbulence occurs, the location of transi- 
tion will greatly affect the character of the wake. A 
detailed discussion of transition to turbulence is 
beyond the scope of this paper. 

Near-wake flow data which are free from support 
interference effects are difficult to obtain. Sting or 
strut mounting of a model in a wind tunnel clearly 
interferes with the base-flow development. Ballistic 
ranges, which have been the premier facilities for 
studying the characteristics of wakes, can provide 
little detail in the base-flow region. Some data 
acquired in magnetic suspension tunnels or with 
wire-supported models exist, but these data are 
generally limited to relatively low temperatures and 
Reynolds numbers. As a result, computational tech- 
niques assume greater importance for providing 
information about this class of flows. On the other 
hand, very little data exist with which to test the 
validity of the computations. 

Computation of high-speed base and near-wake 
flows has received increased attention recently. 
Much of the activity is focused in the area of plumes 
associated with propulsive systems for hypersonic 
flight vehicles. However, nonthrusting base flows 
have received renewed attention. Of particular 
interest are the recent papers of Conti and 
MacCormack, 1 Kim, Loellbach, and Lee 2 and Gnoffo, 
Price, and Braun. 3 In each of these papers, the 
authors numerically solve the Navier-Stokes equa- 
tions in a time-marching manner to obtain a steady- 
state solution for the near-wake flow field. Different 
algorithms are being applied, different gridding 
techniques used, and different gas models employed, 
but the overall solution strategy is the same. The 
same general technique is used in the present paper. 

This work extends the current state of near-wake 
calculations to Reynolds numbers above transition to 
turbulent flow by adding a two-equation turbulence 
model to the equation set. The turbulence mode! 
includes a correction to account for compressibility 
effects. To the authors’ knowledge these are the first 
near-wake computational results including both 
turbulence and nonequilibrium chemistry effects. The 
paper contains several comparisons with perfect gas 
data for base flows and near-wake flows. 


Approach 

The three-dimensional thin-layer Navier-Stokes 
equations are solved using the set of strongly 
coupled, upwind algorithms developed by Molvik and 
Merkle. 4 Both algorithms use upwind differencing, 
Total Variation Diminishing (TVD) techniques in a 
finite-volume framework. The first algorithm is a time- 
marching scheme, and is generally used to obtain 
solutions in the subsonic portions of the flow field. 
The second algorithm is a much less expensive 
space-marching scheme which can be applied in the 
supersonic portion of the flow field. The time- 
marching code has been given the name TUFF (A 
Three-Dimensional, Upwind-Differenced, Finite- 
Volume Flow Solver with Fully Coupled Chemistry), 
and the space-marching algorithm is referred to as 
STUFF. The time-marching scheme uses the zonal 
interfaces described in Ref. 5 to decompose the 
domain into more computationally and geometrically 
efficient blocks. 

The codes currently include perfect gas, 
equilibrium air, and chemical nonequilibrium air 
capability. The equilibrium air model uses the curve 
fits of Srinivasan, Tannehill, and Weilmuenster. 6 The 
species considered in the nonequilibrium air model 
are oxygen (0 2 ), atomic oxygen (O), nitrogen (N 2 ), 
atomic nitrogen (N), nitric oxide (NO), nitric oxide ion 
(NO + ), and electrons (e-). It is assumed that the 
gas mixture possesses a zero net local charge, 
allowing the conservation of electron mass equation 
to be eliminated from the equation set. The reactions 
that are considered are 
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where Mi , M 2 , and M 3 are catalytic third bodies. 
Reaction rates and transport properties are obtained 
from Blottner, Johnson, and Ellis . 7 Wilke’s mixing 
rule 3 is used to compute the mixture viscosity and 
thermal conductivity from those of the individual 
species. Each algorithm is strongly coupled and fully 
implicit, including the chemical source terms. The 
Vigneron, et al.9 approximation is used in the space- 
marching algorithm to allow stable marching in the 
presence of a subsonic viscous wall layer. 

Two turbulence models are included in the code, 
the Baldwin-Lomaxio algebraic model and a two- 
equation k-c model of Nichols, et al.ii. The two- 
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equation model was incorporated for flows with 
multiple turbulent length scales or flows in which 
convection of turbulent quantities is important. The 
low Reynolds number two-equation model is derived 
from the model of Speziale, et a!. 12 The model 
includes compressibility effects on turbulence through 
the compressible dissipation correction of Barker and 
Balakrishnan 13 and through the compressiblity trans- 
formation of Mager 14 in the calculation of y + for wall 
effects. The turbulence equations are solved fully 
coupled with the mean flow equations. 

Results 

Rearward Facing Step 

The code was initially applied to the rearward 
facing step configuration of Smith. 15 The model 
consisted of a 4-in. forward plate and a 12-in. aft 
plate with a 0.443-in. step height. Two cases repre- 
senting the highest and lowest Reynolds numbers 
tested were computed in order to evaluate the code’s 
performance in both laminar and turbulent separated 
flow. Two-dimensional calculations were performed 
for both cases assuming a perfect gas and adiabatic 
walls. A 101 x 71 mesh divided into three zones 
was used for both calculations. The spacing at the 
wall was chosen to correspond to a y + of one. 

The laminar flow case was defined by a Mach 
number of 5.0, a Reynolds number of 1 x 10 5 per 
in., and a total temperature of 680° R. Calculated and 
experimental pressure distributions along the lower 
surface are shown in Fig. 2. Both the base pressure 
and the recompression along the reattachment 
surface are welt predicted by the code. 



10 5 per in., and a total temperature of 620° R. 
Boundary-layer transition location was not measured 
in the test, so the two-equation model was allowed to 
transition naturally. The turbulence model generally 
causes boundary-layer transition to occur prema- 
turely. Based on the short length of the upper plate 
and the relatively low Reynolds number of the 
experiment, the boundary layer on the upper step is 
probably transitional in nature. Calculated and 
experimental pressure distributions along the lower 
surface are shown in Fig. 3. The calculated base 
pressure is slightly greater than the experimental 
result, probably due to the improper treatment of the 
state of the boundary layer on the upper plate, but in 
general the agreement is quite good. 



Fig. 3. Surface static pressure distribution . for a 
turbulent rearward facing step at Mo, = 2.5. 

Laminar Sharp Cone 

Perfect gas calculations of the laminar wake 
behind a sharp, 10-deg half-angle cone in helium at 
Mach 16.35 were performed and compared to the 
experimental data of Murman, et al. 15 and of 
Peterson. 17 The cone had a 1-in. base diameter. The 
tests were performed in the magnetic suspension 
wind tunnel at Princeton and thus are free of support 
interference effects. The Reynolds number was 1.21 
x 10 5 per in. and the total temperature was 560 °R. 
The cone was at 0-deg angle of attack. The solutions 
were found to be sensitive to the viscosity model 
used for the helium gas. The model used for the 
results presented here was 

p = 0.1exp[-1 2.365 + 0.1732in(T) 

- 0.006473ln 2 (T)] (2) 


Fig. 2. Surface static pressure distribution for a Adiabatic no-slip boundary conditions were 

laminar rearward facing step at M» - 5. applied on all walls. The PNS code was used to 

The turbulent flow case was defined by a Mach calculate the flow over most of the forebody. A 91 

number of 2.5, a Reynolds number of 4.6 x x 1 01 computational mesh composed of three zones 
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was used for the TNS code in the near-wake region. 
Calculated density and Mach number distributions 
across the wake are compared to experimental 
results at two stations behind the cone in Figs. 4 and 

5. The calculations are in excellent agreement with 
the data. Features such as the corner expansion and 
the recompression shock can clearly be seen. The 
centerline static pressure distribution is shown in Fig. 

6. The peak recompression pressure is under- 
predicted, causing the pressures downstream of the 
recompression region to be underpredicted. These 
results for the centerline static pressure distribution 
agree well with the calculations of Tassa and Conti. ^ 
The cone base pressure distribution is shown in Fig. 

7. The calculated results agree with the data. The 
base pressure is not constant as is generally the 
case for lower speed base flows. The base pressure 
is low near the corner due to the extremely strong 
expansion of the flow around the corner of the 
model. The base pressure rises towards the model 
centerline because of the recompression which occurs 
when the flow is turned near the wake centerline. 



a. Density 



b. Mach number 

Fig. 4. Distributions across the wake at x/D 
for a sharp cone at Moo = 16.35. 



a. Density 



b. Mach number 

Fig. 5. Distributions across the wake at X/D = 10.0 
for a sharp cone at M„ = 16.35. 



Fig. 6. Centerline static pressure distribution in the 
wake of a sharp cone at M^ = 1 6.35. 


Laminar Blunt Cone 


Calculations of the laminar wake behind a 10- 
percent blunt cone with a 10-deg half-angle were 
performed for the same test conditions as the 
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Fig. 7. Base pressure distribution for a sharp cone 
at Mo. - 16.35. 

previous case and compared to the data of Peter- 
son.^ The blunt cone also had a base diameter of 
one in. The flow over the nose was calculated with 
the time-marching code. This solution was used to 
initialize a space-marching solution over most of the 
forebody. The time-marching code was then run for 
the near-wake region. The space-marching code was 
then initialized in the intermediate wake at an x/D of 
four downstream of the base and marched to an x/D 
of twenty. This was done to demonstrate the space- 
marching capability in the intermediate and far wake. 
Calculated density and Mach number profiles at two 
stations in the wake are compared with experimental 
data at two stations behind the cone in Figs. 8 and 9. 
Both the time-marching (TUFF) and space-marching 
(STUFF) calculated results are shown in Fig. 8. The 
results agree well with the data. Centerline static 
pressure results for both algorithms are compared 
with experimental results in Fig. 10. While calculated 
results overpredict the peak pressure in the recom- 
pression region, there is generally good agreement 
between the calculations and experiment. The cal- 
culations of Tassa and ContP® underpredicted the 


peak recompression pressure for this case. Base 
pressure results are shown in Fig. 1 1 . The peak base 
pressure is overpredicted, but the trends in the experi- 
mental results are reproduced by the calculations. 

Turbulent Sharp Cone 

Perfect gas calculations were performed on the 
turbulent wake behind the sharp, 5.9-deg half-angle 
cone with a rounded corner shown in Fig. 12. The 
flow conditions were for a Mach number of 7.4, a 
Reynolds number of 3.9 x 105 per in., and a total 
temperature of 1,400°R. The computations were 
compared to data obtained in a magnetic suspension 
tunnel at AEDC. The cone was at 0-deg angle of 
attack. The computational grid blocking strategy 
employed with the time-marching code is also shown 
in Fig. 12. 

Boundary-layer transition on the cone was 
assumed to occur following the correlation for sharp 
cones in AEDC tunnels A, B, C, and F 

X tb = 4.6(Re x io-6)-s/7 (3) 

X te = 7.5(Re x i0-6)-5^ (4) 

where X tb is the beginning of transition, X te is the end 
of transition, and Re is the unit Reynolds number per 
inch. Transition was simulated on the forebody with 
the PNS code by damping the eddy viscosity (Pt) 
calculated with the algebraic turbulence model using 
the following: 

p t = 0 ( fi ) 

for x < Xt b and 

p t = (1 -exp[-0.5476(XB - 4.6) 2 ]} ptFT (6) 

for x > X tb , where XB is defined as 

XB = x(Re x 10-5)5/7 (7) 

i S t he fully turbulent eddy viscosity calculated by 
the algebraic turbulence model. The two-equation 
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turbulence model was used after boundary-layer 
transition was completed. The PNS code was used 
to compute the flow over most of the cone forebody. 
A 121 x 91 computational mesh composed of six 
zones (Fig. 12) was used to calculate the flow in the 
near-wake region with the TNS code. The calculated 
and experimental pressure distributions along the 
wake centerline are shown in Fig. 13. The calculations 
accurately predict the cone base pressure and over- 
predict the recompression. Calculated pitot pressure 
distributions across the wake are compared to the 
experimental results in Fig. 14. The major features of 
the computed wake flow field are in excellent agree- 
ment with the measurements. The corner expansion 
and recompression shock are clearly seen. 



a. Density 



b. Mach number 

Fig. 9. Distributions across the wake at X/D = 10.0 
for a blunt cone at Moo = 16.35. 

Blunt Cone with Three Gas Models 

Computational results were obtained for the near 
wake of an 8-deg spherically blunt cone, 127 mm 
long, flying at a velocity of 5,200 m/sec at 0-deg 
angle of attack. Ambient static pressure and tempera- 



Fig. 10. Centerline static pressure distribution in the 
wake of a biunt cone at Moo = 16.35. 
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Fig. 11. Base pressure distribution for a blunt cone at 
Moo = 16.35. 



ture conditions of 100 torr and 300K were chosen as 
typical of AEDC Hypervelocity Range G. A uniform 
cone surface temperature of 2.000K was assumed. 
Computational results were obtained with perfect gas, 
equilibrium air, and chemical nonequilibrium air 
models to assess the effects of chemistry on the 
near-wake flow field. A sketch of the computational 
geometry and solution domain is presented in Fig. 
15. 
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Fig. 13. Wake centerline static pressure distribution 
for a 5.9-deg cone at M*, = 7.4. 

The laminar flow-field solution for the spherical 
nose was obtained with the time-marching code. 
Computation of the flow-field solution over most of 
the afterbody (up to x = 100 mm) was obtained with 
the space-marching code. Instantaneous transition to 
turbulence was assumed to occur at x = 15mm. The 
algebraic eddy viscosity model of Baldwin and 
Lomaxi o was used to approximate the effects of 
turbulence for the first 45 mm after transition, and 
then the k-e model was used for the remainder of the 
calculation over the body and in the wake. 

The near-wake flow field was obtained with the 
time-marching code. The computational domain 
started at x = 100mm and extended 10 base 
diameters downstream from the end of the cone. The 
perfect-gas computation used an 81 x 71 mesh 
divided into two zones. The equilibrium air and 
chemical nonequilibrium air computations used a 121 
x 101 mesh composed of two zones. The two zones 
used for the afterbody and near-wake computation 
are shown in Fig. 15. Constant temperature, no-slip 
wall boundary conditions were specified along the 
base. For the nonequilibrium computation, a non- 
catalytic wall boundary condition was used on all 
body surfaces. Local time-stepping was used for all 
of the cases considered. 

Pressure contours, normalized by the free-stream 
static pressure, and velocity vectors for the perfect- 
gas case are presented in Figs. 16 and 17, 
respectively. The principal features of the near-wake 
flow are apparent. The corner expansion, recom- 
pression shock, primary recirculation region, and 
viscous boundary layer on the base all are clearly 
discernible. In addition, secondary and tertiary recir- 
culation regions can be seen near the corner of the 
cone. An embedded shock can be seen near the 





c. X/Dbase - 5- 02 

Fig. 14. Pitot pressure distributions across the wake 
for a 5.9-degree cone at = 7.4. 

body axis and close to the base. The primary recircu- 
lation region is locally supersonic near the body axis 
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y ob 

Fig. 15. Eight degree spherically blunt cone. 



Fig. 16. Pressure contours (P'P<») for a blunt cone with a 
perfect gas model. 



Fig. 17. Velocity vectors for a blunt cone with a perfect gas 
model. 


and is compressed through the embedded 
shock from the presence of the base. Conti 
and MacCormacki and Kovenya and Lebedevis 
reported similar features. Comparisons of the 
wake centerline pressure and temperature 
distributions obtained with the different gas 
models are presented in Figs. 18 and 19. The x 
coordinate is measured from the base of the 
cone. The wake stagnation point occurs 
approximately 1.0 to 1 .5 base diameters 
downstream from the base of the cone 
depending on the gas model used. Downstream 
of the wake stagnation point the flow expands 
to a pressure below the free-stream pressure. 
Equilibration of the pressure with the free 
stream will occur over a downstream distance 
of many base diameters. Compression of the 
flow through the embedded shock also 
generates centerline temperatures that are well 
above the specified cone surface value. 
Additionally, the choice of the gas model greatly 
influences the peak centerline temperature. The 
perfect-gas model predicts a peak temperature 
of approximately 19 times the free-stream 
value, the nonequilibrium model predicts a peak 
temperature that is approximately 13.5 times 
the free-stream value, and the equilibrium 
model predicts a peak temperature that is 
approximately 20 percent less than the 
nonequilibrium value. These trends are 
expected because of the targe amounts of 
energy absorbed by the chemical reactions. 
After 10 base diameters, the temperature is 10 
to 12 times greater than the free-stream value. 
The wake centerline electron density distribu- 
tion, an important parameter to the radar 
signature, is presented in Fig. 20. The peak 
electron density occurs in the hot core of the 
wake and decreases as the flow accelerates 
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and cools downstream of the wake stagnation point. 
The calculated effects of the various gas models 
agree with the trends reported by Kim, Loellbach, 
and Lee. 2 




Fig. 19. Wake centerline temperature distributions for 
a blunt cone with three gas models. 


Summary 

A computational technique applicable to high- 
speed base flows has been demonstrated. The 
technique is based on solving the three-dimensional 
thin-layer Navier-Stokes equations for either a perfect 
gas, equilibrium air, or a chemical nonequilibrium air 
model. An efficient space-marching algorithm is used 
whenever possible on the body and in the far wake. 
A blocked time-marching algorithm is used to 
calculate the embedded subsonic flow in the base 
region. A two-equation turbulence model which 
includes compressibility effects is used to extend the 
approach to high Reynolds numbers. 


A series of computations has been made to 
compare with base-flow and near-wake flow data. 
These test cases have free-stream Mach numbers 
ranging from 2.5 to 16.35 and encompass both 
laminar and turbulent flow conditions. The experi- 
mental data used for comparison were obtained in 
wind tunnels with stagnation temperatures low 
enough that no dissociation of the test gas was 
present. Good agreement with the data was obtained 
in ail cases. 

Further demonstration of the computational 
capability has been accomplished by carrying out 
calculations for a blunt cone configuration with free- 
stream conditions typical of those encountered in a 
ballistic range. The Reynolds number was assumed 
to be high enough to cause boundary-layer transition 
on the cone frustrum. These computations were 
performed assuming a perfect gas, equilibrium air 
and chemical nonequilibrium air. 
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ABSTRACT 

This paper presents the results of a feasibility 
study of a hydrocarbon scramjet design utilizing 
an augmented preburner upstream of the main 
fuel injector locations. The combustor design 
evaluated here is for a small hypersonic research 
vehicle. It consists of a preburner into which a 
small amount of fuel is burned with on-board 
liquid oxygen and injected into the main airflow, 
upstream of the main fuel injector locations, thus 
ensuring that combustion is present and 
uninterrupted. Two degrees of analysis are 
presented including a one-dimensional cycle 
analysis and a complete computational fluid 
dynamic analysis with finite-rate chemistry and a 
two-equation turbulence model. Comparison of 
these analyses show good agreement when the 
CFD-predicted fuel consumption schedule is 
used in the cycle analysis. 

INTRODUCTION 

The advent of the National Aero-Space Plane 
(NASP) has generated renewed interest in 
hypersonics research. The U. S. Air Force and 
NASA are pursuing the design and development 
of a single-stage-to-orbit (SSTO) vehicle in the 
NASP X-30. Both Europe and Japan have also 
proposed hypersonic vehicle design activities. 
Applications of hypersonic vehicle concepts to 
high-speed missions, both military and civilian, 
are also underway in government and industry. 
There exists an increasing need for hypersonic 
research vehicles (HRV) to demonstrate 
integrated aerodynamic, propulsion, and 
structural technologies for hypervelocity design 
and to develop a research database for reducing 
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the risk involved in the development of 
operational hypersonic vehicles. Conceptual 
design activities are currently under way at NASA 
Ames Research Center to determine the 
feasibility of such a research vehicle utilizing 
near-term technology (Fig. 1). 

The research goals of this activity are to provide 
an understanding of the underlying physics, 
verification of design tools, and validation of the 
technologies and systems needed. The 
research requirements are classified into two 
areas: 1) basic research, and 2) systems 
technology demonstration, which addresses 
programmatic research issues and overall vehicle 
system integration and performance. The 
disciplinary research requirements include 
aerodynamics and aero-thermodynamics of 
hypersonic flight, hypersonic air-breathing 
propulsion system performance, structures and 
materials characteristics, and finally 
instrumentation requirements for hypersonic 
vehicles. 

The propulsion system for this HRV study is a 
hydrocarbon supersonic combustion ramjet 
(scramjet). Both NASA and DoD are currently 
reviewing a wide range of missions requiring a 
high-speed performance capability [1-3]. For a 
Mach number range of 4 to 10, hydrocarbon 
fuels provide sufficient engine specific impulse 
(Isp) performance, heat sink capability, and offer 
the potential of reduced vehicle size compared 
with hydrogen-powered designs. In addition, 
the handling and infrastructure requirements for 
the hydrocarbon fuels have a distinct advantage 
compared to cryogenic hydrogen. 

The slow reaction rates of a hydrocarbon/air 
mixture in a supersonic stream can have a 
significant impact on the inlet/combustor design. 
Because of the size limitation of an engine on a 
small research vehicle, a mechanism is required 
to provide sufficient fuel/air temperatures for 
burning within the combustor. The concept of 
employing an in-stream, embedded ramjet as a 
pilot light has been proposed [3-4] and appears 
to be a promising technique for maintaining 
combustion. An alternative is to use a liquid 
oxygen (LOX) augmented pre-burner located 
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upstream of the main fuel injectors to promote 
burning in the combustor. Because the LOX 
would be stored onboard the research vehicle, 
the required pre-burning should be kept to a 
minimum to reduce the impact on the vehicle 
design and gross weight. For the relatively short 
hypersonic research mission cruise times (5 to 1 0 
minutes), the impact of this additional onboard 
mass will not have a significant impact on the 
HRV design. 

In this paper, the feasibility of a hydrocarbon 
scramjet engine utilizing a LOX/preburner 
concept is addressed. An engine design for a 
small hypersonic vehicle is proposed and then 
analyzed. Two degrees of analysis are 
presented in this paper. A one-dimensional 
cycle code is first used to define a baseline 
configuration and the operating conditions of the 
engine. A computational fluid dynamic (CFD) 
analysis is then performed to predict the details 
of the engine flow field. The results of the CFD 
analysis are then fed back into the cycle code for 
the purpose of further refining the design and 
computing an overall vehicle performance. 
Comparisons between the cycle and CFD results 
are also included. 

CYCLE ANALYSIS 

As opposed to accelerator-type missions (e.g. 
SSTO) where the mass capture characteristics of 
the vehicle are most important, cruise 
configurations place an emphasis on the 
aerodynamic performance of the design. Hence 
vehicle lift-to-drag ratio (L/D) or the product of 
L/D and engine Isp (a parameter proportional to 
Brequet range factor) becomes more important. 
A waverider configuration, with high hypersonic 
L/D potential, was selected as the baseline 
configuration for the present study of the HRV. 
The configuration/engine installation was 
numerically optimized to maximize (L/D) x Isp, 
using forebody shape, ramp angles and cowl 
position as optimization parameters. For this 
analysis, the forebody is a Mach 8, waverider 
configuration and the ramp angles and cowl 
position were designed to produce shock-on-lip 
and shock-on-shoulder[5]. The engine 
geometry and operating parameters were held 
fixed in the optimization procedure. The HRV 
conceptual design analysis and sizing was 
computed using the hypersonic vehicle 
synthesis code (HAVOC) of Ref. 6. The 
numerical optimization was performed using the 
methods of Ref. 7. 


As part of the HRV vehicle synthesis, the nose- 
to-tail propulsion module of the HAVOC code 
was used to predict installed engine performance 
of the hydrocarbon scramjet. The forebody and 
nozzle flow fields are computed using an inviscid 
2-D real gas, weak wave/oblique shock analysis. 
In addition, the nozzle flow field was computed 
assuming frozen flow (mole fractions taken at the 
combustor exit plane). The combustor flow is 
solved using the one-dimensional mass, 
momentum, and energy equations for the 
fuel/air mixture, and marching through the 
combustor with a specified number of steps. 
Multiple fuel injection stations and LOX- 
augmented preburning injection are accounted 
for. Overall engine heat balance is computed 
using an input combustor skin friction coefficient 
as a function of freestream Mach number. For 
the initial analysis, the combustion efficiency was 
taken from Ref. 8, which accounts only for mixing 
efficiency, with no reaction delay time effects 
included. Modeling the combustion efficiency as 
a function of location in the combustor (i.e. a heat 
release schedule) enables the one-dimensional 
code to predict flow properties and engine 
performance using a multi-step approach. The 
combustor efficiency model must come from 
either experimental data or more detailed 
calculations. For results presented later, an 
improved combustor efficiency was computed 
with CFD and then implemented in the cycle 
analysis code. 

CFD ANALYSIS 

The design of future high-speed propulsion 
systems such as those for NASP and the HRV 
will for the most part be based on CFD. The 
need for CFD in scramjet propulsion design 
stems from the fact that scramjet propulsion is 
difficult to test in ground-based facilities and has 
yet to be demonstrated in flight. Therefore, CFD 
will continue to play a major role in the design of 
scramjet propulsion systems. 

The requirements of a flow solver for hypersonic 
scramjet computations are very demanding. The 
flow solver must be capable of predicting the 
three-dimensional flow of a highly-turbulent 
mixture of reacting gasses with separated 
regions and strong flow field discontinuities. 
Upwind-differenced schemes offer an appealing 
approach to solutions of scramjet flow fields 
because of their ability to capture strong flow 
field discontinuities without user-specified 
smoothing terms. Further, multi-equation 
turbulence models become a requirement for 
accurately computing the turbulent shear layers 
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and boundary layers in a scramjet engine. Finally, 
a strongly coupled, nonequilibrium chemistry 
capability is required to compute the highly- 
reactive combustion processes that are present 
in a hydrocarbon scramjet engine. 

Numerous unsteady numerical methods have 
recently been developed [9-18] that address 
chemical nonequilibrium effects in both internal 
and external flows. The numerical scheme 
chosen for this study is the time-marching flow 
solver, TUFF, originally developed by Molvik and 
Merkle[9] for external hypersonic, reacting flow 
fields. Even though the TUFF code was 
originally developed for external hypersonic 
computations, it has many of the characteristics 
required for internal, scramjet computations. It 
employs a finite-volume philosophy to ensure 
that the scheme is fully conservative. The 
inviscid fluxes are obtained by employing a 
temporal Riemann solver that fully accounts for 
the presence of a multicomponent mixture of 
gasses. A total-variation-diminishing (TVD) 
technique of the type outlined by Chakravarthy 
[19] allows extension to higher orders of 
accuracy without introducing spurious 
oscillations. The scheme employs strong 
coupling between the fluid-dynamic and 
chemistry equations. It solves the thin-layer, 
Reynolds-averaged, Navier-Stokes equations 
with viscous terms in all directions and has 
generalized boundary conditions. The scheme 
is made implicit with full linearization, approximate 
factorization and by employing a modified 
Newton iteration process to reduce linearization 
and approximate factorization errors. The 
turbulence model included in the original code is 
the algebraic model of Baldwin and Lomax [20]. 
Recently, a zonal scheme was incorporated [21] 
that utilizes patched interfaces to maintain 
conservation and spatial accuracy. For a 
complete description of the governing equations 
and numerical scheme, the reader is referred to 
Ref. 9. 

Several modifications were required to the 
original TUFF code for scramjet computations. It 
was necessary to change the kinetics mode! from 
pure air to hydrocarbon / air including the reaction 
rates and mass action coefficients. Similar 
changes to the thermodynamic and transport 
data were also required. Further, a higher-order 
turbulence model with compressibility effects 
was added to improve the prediction of the many 
viscous phenomena. Finally, a boundary 
condition procedure was added to model the 
injection ports. These modifications are detailed 
below. 


Kine ti cs Mo del 

Flow field analysis of complete finite-rate 
chemical kinetics for hydrocarbon based fuels 
with air can be very complex and resource 
exhaustive. Hydrocarbon/air models have been 
proposed with tens of species and hundreds of 
reactions. CFD codes that include these 
complete models are currently limited to 1-D 
simply because of the exhaustive computer 
resource requirements for extension to higher 
dimensions. Simplified reaction mechanisms for 
the hydrocarbon combustion process offer a 
solution to this problem by foregoing the details 
of the combustion process. In this approach, 
various species and reactions are combined and 
simplified while preserving the net effect of the 
reaction processes. Since kinetic details are not 
required in the present study, a simplified 
reaction mechanism was sufficient. 

For the scramjet propulsion system analysis, the 
two-step reaction mechanism of Westbrook and 
Dryer [22] was employed to simulate the 
combustion of fuel with air. For hydrocarbon 
scramjet propulsion, the liquid fuel can be used 
as a coolant on various portions of the aircraft. 
This results in an endothermic reaction that can 
vaporize and dissociate the liquid fuel. Gaseous 
ethylene was used in this analysis, as a surrogate 
fuel intended to represent the products of this 
endothermic reaction. The kinetics model for 
ethylene is outlined below. 


C 2 H 4 + 20 2 <-> 2CO + 2H 2 0 (1 ) 

CO +2 0 2 h C0 2 (2) 

The first reaction represents the combustion of 
ethylene with oxygen. The second reaction is 
included to complete the combustion process 
and to serve as an equilibrium mechanism for the 
products that then improves the predicted heat 
of reaction and adiabatic flame temperature. This 
equilibrium step is especially necessary since 
substantial amounts of CO and H 2 can exist in 
the combustion products along with CO 2 and 

h 2 o. 

The values of the forward and backward reaction 
rate constants are evaluated with Arrhenius fits to 
existing laminar flame data. These take the form: 

k = AT n exp(-E a /RT) (3) 


3 



where the pre-exponential factor, A, and the 
activation energy, E a , are constants for a 

particular reaction. The temperature- 
dependence exponent, n, was set to zero for ail 
of the reaction rate fits. The reaction rate 
constants recommended by Ref. 22 are 
tabulated below: 


Table 1 Reaction Rates 


Eauation 

A 

Ea 

If 

2.4 x10 12 

30 

1b 

0.0 

- 

2f 

4.0 x10 14 

40 

2b 

5.0x10 8 

40 


The backward reaction rate for equation 1 was set 
to zero. The units for the pre-exponential 
constant and for the activation energy are cm- 
sec-mole-kcal-kelvins. 

The mass production rate of each species is 
determined with a modified law of mass action in 
this simplified model. The reaction rates 
evaluated with the modified law, are shown 
below: 

R, =[C 2 H 4 J 0 ' 1 [O 2 ] 1 ' 65 k 1f (4) 

R 2 - [C0]1'°[H 2 0]° 5 [0 2 ] k 2f - [C0 2 ] 10 k 2b (5) 

with the bracketed terms representing the molar 
concentrations of each indicated species. The 
mass production rate for each species, s, can 
now be determined with the following: 

cb s =^s[(v1s-v!s)Rl+(v^s-v^s)R2] 

where, vfcs and v £s are the reactant and product 
stoichiometric coefficients and M s is the 
molecular weight of species s. Throughout this 
kinetics model, molecular nitrogen, N 2 , has been 
treated as an inert species. 

Thermodvnamic/Transport Properties 

The thermodynamic and transport properties for 
the individual species in the hydrocarbon/air 
mixture were obtained from Ref. 23. For these 
curve fits, the properties are expressed solely as 
functions of temperature. The temperature 
range for this data is typically limited to 300- 
5000°K. However, research is currently 


underway at NASA Lewis Research Center to 
expand this range [24]. For details on the 
computation of mixture properties in the TUFF 
code, see Ref. 9. 

Turb ule nc e M odeli ng 

Modeling scramjet flow fields with CFD requires 
an advanced turbulence model capable of 
accurately accounting for compressible turbulent 
shear layers and jets. This was accomplished in 
the present study with the incorporation of the 
low Reynolds number K-e turbulence model 
originally developed by Jones and Launder [25]. 
The compressibility correction of Zeman [26] was 
also included in the two-equation formulation to 
improve the computation of compressible shear 
layers. The modified numerical constants were 
suggested by Ref. 26 to produce better 
boundary layer results as well as improve shear 
layer growth rates. The entire turbulence model 
was transformed to a generalized finite-volume 
coordinate system and strongly coupled with the 
existing flow solver, including the viscous and 
inviscid flux computation and the source term 
treatment. 

Injector Boundary Condition Procedure 

For the engines studied in this paper, all of the 
injectors {including preburner and main) were 
designed to be supersonic in the boundary- 
normal direction. This simplifies the boundary 
condition procedure since merely specifying of 
the injector variables is required for a supersonic 
inflow boundary condition. This type of injector is 
also quite practical in a scramjet since injector 
pressures are typically high enough to choke the 
injector flow. The supersonic inflow boundary 
condition was imposed only on those cell faces 
that correspond to an injector exit. No-slip, 
viscous boundary conditions were imposed on 
the cell faces adjacent to the injector exit. This 
led to the use of rectangular injectors to avoid 
further complication of the boundary condition 
procedure and the grid generation process. 


RESULTS 

Three sets of results are included in this paper: a 
shear layer test case, and two- and three- 
dimensional scramjet results. Because of the 
absence of high-speed, hydrocarbon 
combustion data, a validation case for the 
combustion model of Ref. 22, under scramjet 
conditions, is not possible. Validation for this 
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model Is therefore limited to the laminar flame 
comparisons with shock tunnel data by the 
originators of the model. High quality, CFD 
validation experiments in the area of high-speed 
hydrocarbon combustion are therefore greatly 
needed. 

Sh ea r. L ayer Js s L Casa 

The first test case is that of two parallel streams of 
perfect gasses that are initially separated by a 
splitter plate with boundary layers of zero 
thickness (Fig. 2). This test case is intended to 
validate the Jones and Launder turbulence 
model, with the Zeman compressibility 
correction, for high-speed shear layers. The 
flow conditions were chosen to correspond with 
the computation of Viegas [27] in which a 
comparison with experiment was conducted. In 
this test case, the total temperature and static 
pressure of the two air streams were matched 
and set to 833°K and 1 atmosphere, 
respectively. The Mach numbers of the two 
streams were 4.92 and 1.1. Strong 
compressibility effects are expected for these 
flow conditions. 

The results of the shear layer test case are shown 
in Figs. 3-5. Note that for the figures presented 
in this section, the high Mach number stream is 
on the top (positive y direction). Results using 
the standard Jones and Launder K-e turbulence 
model with no compressibility correction are 
presented first. A direct comparison of results 
with those of Viegas are shown in Figs. 3 and 4. 
These figures present profiles of velocity and 
turbulent kinetic energy profiles at five axial 
stations in the shear layer, plotted as a function of 
distance across the shear layer normalized by the 
local vorticity thickness of the shear layer. These 
results indicate a linear growth of the shear layer 
with axial distance and also suggest that a profile 
similarity exists for a shear layer with zero 
boundary layer thickness at the end of the splitter 
plate. These figures also show that the shear 
layer penetration is much greater into the low 
Mach number stream. There is also a general 
turning of the high-speed stream into the low- 
speed stream. This is in agreement with the 
results of Viegas. 

Comparison of these results with those of Viegas 
are excellent except for a slight vertical drift that 
can be seen in both Figs. 3 and 4. This 
difference is attributed to the far-field boundary 
conditions used by Viegas. Those boundary 
conditions produced a small vertical velocity 
throughout the flow field that tended to further 


turn the shear layer into the low Mach number 
stream. Comparisons of the profiles are 
otherwise in direct agreement including the peak 
turbulent kinetic energy, the velocity profile and 
shear layer thickness. 


The linearity of the growth rate is easily seen in 
Fig. 5. This figure shows the vorticity thickness 
as a function of the axial distance from the splitter 
plate. Vorticity thickness is defined as follows: 


8 vv = (Ui-U 2 )/| 


I'au'l 

3y 


V yjJ 7m ax ( 7 ) 

Comparison of the standard Jones and Launder 
two-equation turbulence model with that of 
Viegas shows very good agreement. The 
predicted growth rate for the two-equation model 
with the compressibility correction is also plotted 
here. This plot shows a reduction in the growth 
rate from the standard model which has been 
shown by Refs. 26 and 27 to be in better 
agreement with experiment. 


Two-Dimensional Scramjet Result 


As a first attempt to determine the feasibility of a 
hydrocarbon scramjet engine with augmented 
preburning, a two-dimensional analysis was 
conducted. The initial geometric definition for 
the scramjet engine, including throat height, 
shock isolator length, and combustor length, 
and combustor area ratio were taken from Ref. 4. 
The embedded ramjet section was removed and 
a hydrocarbon/LOX preburner was added. A 
mixing section aft of the preburner station was 
also added to allow mixing of the preburner 
exhaust gasses with the oncoming air so as not 
to suffocate the main burner jets of oxygen. A 
schematic showing the scramjet concept is 
presented in Fig. 6. Backward facing steps, prior 
to the main fuel injection station were included to 
act as flame holders and mixing enhancement 
mechanisms. Finally, for the purpose of 
comparison, a preburner was only employed on 
the top of the engine leaving the lower main 
burners exposed to just the oncoming air stream. 

Preliminary HAVOC design results at a Mach 
number of 8 and a dynamic pressure of 1500psf 
for the waverider HRV with two ramps and with 
both shock on shoulder and on cowl lip indicated 
that a contraction ratio of roughly 14 was 
achievable. Details of this inlet design are 
included in Ref. 5. Using a guideline of 
approximately 1000°K as auto-ignition of a 
gaseous ethylene/oxygen mixture, the 1 -D cycle 
code was run with the LOX augmentation 
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preburning option to compute parametrically, the 
required fuel and oxygen flow to achieve an 
equivalent mixed 1-D temperature at the main 
fuel injector station equal to the auto-ignition 
value(Fig. 7). For an engine with an equivalence 
ratio of 1 (stoichiometric) this resulted in roughly 
2.5% of the fuel being directed to the preburner 
which was then burned stoichiometricaly with 
onboard LOX. A heat balance on the vehicle and 
engine was used to compute the fuel total 
temperature. The preburner and main burner 
pressure and velocity were then selected to 
produce supersonic normal injection and the 
required exit areas were also computed. The 
resulting operating conditions of this engine are 
given in Table 2. 


Table 2 . Engine Operating Conditions 



Inlet 

Preburner 

Main 

Gas 

Air 

Products 

C 2 H 4 

Mach No. 

3.84 

2.5 

2.5 

T(°K) 

795 

2897 

403 

Ptetm) 

.95 

1.76 

1.73 

Areafin) 

2.0000 

0.0268 

0.0820 

Anale 

- 

90° 

90° 


The area in the above table is based on a one 
inch width section. 


The CFD results for this two-dimensional 
approximation are shown in Figs. 8-12. 
Throughout the CFD simulation process, the 
ingestion of a thick boundary layer formed on the 
forebody of the hypersonic vehicle by the 
scramjet was neglected. The grid for this 
computation spanned 150 cells in the axial 
direction and 74 cells in the cross flow direction 
and was generated algebraically. For this 
turbulent computation, the grid spacing at the 
wall was set to 2.5 x I0' 5 m. The pressure 
contours for the entire scramjet engine are 
shown in Fig. 8 and indicate a somewhat smooth 
pressure variation throughout the engine except 
in the vicinity of the injectors and steps. The 
preburner exhaust gasses produce a shock that 
traverses the entire height of the scramjet. The 
pressure rise of this shock is then quickly 
suppressed by a weak expansion that is caused 
by a slight geometric expansion region starting 
just after the preburner location. These 
preburner gasses then are convected 
downstream and seem to adhere to the upper 
surface for this two-dimensional computation. 

The upper and lower step flows exhibit entirely 
different behavior, because the upper step was 


exposed to the hot preburner exhaust gasses 
and the lower step was not. The lower step flow 
field exhibits a steady behavior and contains 
many of the features that are typical of a reward 
facing step flow field. A recirculating region 
located aft of the step is present but is slightly 
enhanced by the existence of the normal fuel 
injector placed one quarter of a step height 
downstream of the step. The gas behind this 
step is entirely comprised of fuel. A 
recompression shock is also clearly visible in the 
pressure contours of Fig. 8. Combustion of the 
lower gasses begins about five step heights aft 
of the step and is visibly enhanced by the 
recompression shock. This enhancement effect 
can be seen about 8 step heights aft of the step 
in the temperature contours of Figs. 9 and 10. 

The upper step on the other hand produces a 
violently unsteady flow field. The hot preburner 
exhaust gasses on the upper surface interact 
with the air and fuel to produce intermittent 
periods of combustion aft of the step. The 
predicted flow field is highly irregular and no 
cyclic behavior was observed. It should be noted 
that a local time-stepping routine was used and a 
time-accurate scheme could produce an entirely 
different solution. At the point in the 
computation that these figures were generated, 
a pure fuel jet was observed that penetrated well 
into the oncoming air stream. This phenomenon 
can be seen both in the fuel contours of Fig. 1 1 
and the velocity vectors of Fig. 12. 

Several deficiencies emerged in this two- 
dimensional analysis of a clearly three- 
dimensional problem. First, the geometry 
required modification to accommodate a two- 
dimensional computation. This consisted of 
considerably shortening the injection port 
lengths to maintain a constant injection area as 
the ports were changed from series of individual 
ports to a single slot. Second, the jet penetration 
and mixing efficiency are significantly reduced in 
a two-dimensional computation. This is the most 
prevalent effect and can significantly affect the 
predicted performance of the scramjet engine. 
And finally, any side wall effects are entirely 
neglected in a 2-D computation. For these 
reasons, a three-dimensional analysis was 
performed. The three-dimensional CFD analysis 
presented in the following discussion addresses 
the first two issues stated here, although it still 
neglects any side wall effects. 
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Three-Dimensional Scramiet Result 


Several modifications were made to the two- 
dimensional scramjet engine design before the 
three-dimensional computation was performed. 
First, the forebody geometry was further 
optimized with the HAVOC code, resulting in 
slightly different inlet conditions. These new 
inlet conditions are included in Table 3. Next, 
the two-dimensional results indicated that the 
steps filled up entirely with fuel and did not 
serve the purpose of a flame holder. Further, 
they produced a negligible effect on mixing, 
and cycle analysis indicated that engine 
performance is better served with this aft-facing 
area distributed in the main combustion region. 
For these reasons, the steps were eliminated in 
the three-dimensional scramjet design. The 
flow of LOX was then reduced to the preburner 
resulting in a fuel-rich preburner. This reduced 
the required amount of onboard LOX with only a 
minimal effect on temperature distribution prior 
to main fuel injection. This is because the fuel- 
rich preburner exhaust gasses continue to react 
with the air after injection. Finally, preburning 
was employed both on the top and bottom of 
the scramjet to improve engine efficiency. Fig. 
13 shows a schematic of the modified scramjet 
engine. 


As in the two-dimensional design, the preburner 
and main burner injectors were designed to 
produce supersonic injection. A further 
consideration in the three-dimensional design 
was the penetration distance of the injectors. 
Reference 28 gives a relationship for the jet 
penetration distance as a function of the jet and 
freestream Mach numbers, the momentum ratio 
of the two streams, the angle of injection and the 
jet diameter, dj. This relationship is given below: 
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The preburner injection was designed to only 
penetrate through the boundary layer , whereas 
the main injection was designed to reach well into 
the air stream for better mixing. This resulted in a 
guideline of h=0.5cm for the preburner and 
h=2.0cm on the main burner. Because of the 
iarge amounts of fuel being injected through the 
main injectors at stoichiometric conditions, the 
main fuel injectors were designed with a 
streamwise length-to-width aspect ratio of 5 and 


an angle of injection of 30 degrees to help 
reduce blockage. The injectors were laterally 
spaced one inch apart on both the top and 
bottom of the scramjet. The top injectors were 
then offset one-half inch to a produce a 
staggered injection for the purpose of increasing 
jet penetration and to avoid the additional losses 
of impacting jets. The preburners were aligned 
with the main fuel injectors to ensure that the hot 
preburner gasses fell in near vicinity of the main 
fuel. The HAVOC code was run with this new 
engine and produced the following engine 
operating conditions. 

Table 3 . Engine Operating Conditions 




Preburner 


Gas 

Air 

Products 
and C2H4 

C 2 H 4 

Mach No. 


1.1 

2.2 

■QBSH 

833 


801 

Pfatm) 

.86 


2.49 

Areattn) 

Em 

■•mi 


Anale 

- 

90° 



A three-dimensional CFD analysis of this engine 
was performed. Because of the periodicity of 
the engine in the absence of any side walls, only 
the flow between the centerline of the top 
injectors and the centerline of the bottom 
injectors was actually solved. The grid for this 
computation was generated algebraically and 
contained 119 cells in the axial direction, 60 cells 
from top to bottom, and 16 cells laterally. The 
grid spacing on the walls was set to 4.0 x 10‘5m. 
All of the injector exits were modeled as 
rectangles containing 8 cells in each of the axial 
and lateral directions. This three-dimensional 
computation required 1850 iterations leading to 
1 02 hours on a Cray-YMP. The computation was 
halted after no plotable difference was seen with 
further iteration. 

The results of the three-dimensional CFD 
analysis are shown in Figs. 14-26. The pressure 
contours on both the symmetry plane containing 
the top injector and the one containing the lower 
injector are shown in Fig. 14. As in the two- 
dimensional results, the pressures are smoothly 
varying except in the vicinity of the injectors. 
The large influence of these injectors is clearly 
visible in these figures and affects the entire flow 
path. Shocks emanating from both the 
preburner and main fuel injectors traverse the 
height and width of the computational space. 
These shocks can significantly affect the 
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efficiency of the engine, and any further 
refinement of this design will address the losses 
associated with these large structures. The 
temperature contours (Fig. 15) show the injector 
penetration and the temperature rise caused by 
combustion. This figure indicates that the 
penetration of the combustion region is slightly 
more than half the height of the scramjet. A 
significant degree of penetration is present, yet 
avoids adverse effects such as jet-jet 
interactions and jet-wall interactions. 

The depth of penetration and degree of 
combustion are better seen in the fuel and 
oxygen mass fraction contours of Figs. 16 and 
1 7. The fuel contours clearly show the regions 
of unburned gaseous fuel and the oxygen 
contours show the locations that are deficient in 
oxygen indicating the penetration distance of 
fuel and products. These oxygen contours 
show that the injector gasses do indeed 
penetrate far into the flow path. Also shown in 
the fuel contours are the fuel-rich preburner 
gasses The depletion of the fuel and oxygen in 
the preburner gasses indicates that combustion 
is present and exhausts nearly ail of the fuel in 
the preburner gasses. 

The velocity vectors of Figs. 18 and 19 show 
some details of the injector flow fields. The 
velocity vectors in the immediate vicinity of the 
preburner exit exhibit a spreading behavior that 
is typical of an under-expanded jet. This 
phenomenon is also present in the main injector 
region but is less visible because of the 
inclination of the vectors. A separation region is 
present in the preburner injector region that 
reaches seven jet diameters upstream of the 
injector. This phenomenon is absent near the 
main fuel injectors because of the reduced 
angle of injection. Also shown in this figure is 
the increased penetration that can be achieved 
with normal injection. 

Figures 20,21 and 22 contain crossflow contour 
plots of temperature, fuel and water 
respectively, at various axial locations. The axial 
locations correspond to the following: 1) the 
back of the preburner station, 2) just upstream of 
the main fuel injection, 3) the back of the main 
injector station, 4) within the combustion 
chamber, and 5) the combustor exit. The exact 
axial locations are indicated on the plots. 

The temperature contours of Fig. 20 clearly 
show the mechanism that is studied in this 
paper. Fig. 20(a) shows the hot preburner 
gasses that emerge from the preburner injector 


ports. These gasses mix and react with the 
oncoming air stream but still contain a very hot 
core just before main fuel injection (Fig. 20(b)). 
This hot core, falling just above the main fuel 
injection, serves as a "pilot light" for main fuel 
injectors causing combustion of the main fuel to 
instantaneously occur (Fig. 20(c)). The main 
fuel injectors were designed to produce a 
significant amount of penetration without 
traversing the entire height of the scramjet. This 
was accomplished and is clearly shown in the 
combustion chamber temperature contours 
(Figs. 20(d&e)). These figures indicate that the 
concept of preburning does indeed accomplish 
the task of maintaining combustion at the main 
fuel injection station and that an injector can be 
designed to provide significant flow path 
penetration without unstarting the engine. 

The effect of the upper surface corner on 
combustion is shown in Figs. 21 and 22. An 
expansion wave emanating from this corner 
causes the density and temperature to decrease 
having an adverse effect on the rate of 
combustion and mixing. This wave affects the 
upper gasses before reaching the lower gasses. 
Therefore the expansion has a greater effect on 
the upper gasses. For this reason, there are 
more unburned and unmixed gasses present in 
the upper region of the scramjet. 

The ability of the 3-D CFD finite-rate model to 
provide the equivalent 1-D combustor efficiency 
(i.e. the heat release schedule) allows the 1-D 
cycle analysis code to predict thermodynamic 
flow properties through the combustor, and 
hence compute engine and vehicle 
performance. Sensitivity analysis using the 
cycle code indicates that at Mach 8 a 1% change 
in overall combustion efficiency represents 
approximately 1% change in cowl-to-tail Isp and 
0.8% change in axial thrust coefficient. Hence 
having a good combustor efficiency model 
permits the use of a 1-D cycle code to model 
engine performance with sufficient accuracy for 
conceptual design analysis. 

The remaining figures present a comparison of 
the CFD results with those of a 1-D cycle 
analysis. The combustor efficiency computed 
by the CFD solution was implemented in the 
cycle code since no other schedule was 
available for this engine design. This was 
accomplished by curve fitting the average fuel 
fraction schedule (Fig. 23) and using this 
schedule in the 1-D cycle code. Both the CFD 
predicted schedule and the curve fit are shown 
in Fig. 23. Also shown on Fig. 23 is the 
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predicted amount of carbon monoxide from both 
the CFD and 1-D cycle analysis. The CFD 
analysis predicts a higher amount than the 1-D 
cycle code. This is caused by the difference in 
the equilibrium mechanisms of the two codes 
and suggests an improvement to the simplified 
kinetics employed in the CFD solver for scramjet 
computations. Comparison of the average 
temperature, the momentum-averaged pressure 
and the mass-averaged velocity from the CFD 
analysis with the results of the cycle code show 
general agreement (Figs. 24-26). The 
discrepancies can be attributed to the improved 
ability of the CFD method to account for detail 
and the difference in the equilibrium 
mechanisms. 

CONCLUSIONS 

As a result of this effort, two methods for the 
conceptual design and analysis of hydrocarbon 
fueled scramjets were developed and 
demonstrated: a cycle code with a simplified 
nose-to-tail flow field analysis capability, and a 
complete 3-D CFD code with finite-rate chemistry 
for hydrocarbon/air combustion. Utilizing the 
heat release schedule predicted by the CFD 
analysis in the cycle code, the combustor flow 
field properties and level of combustion product 
constituents computed by the two respective 
methods agree. Verification of the 1-D cycle 
code results by the CFD analysis encourages 
further application of the cycle code to a broader 
range of engine design parameters. Promising 
configurations can then be analyzed in detail 
using CFD. However, there is currently a void of 
high-quality, CFD validation-type data for high- 
speed hydrocarbon combustion, leaving an 
uncertainty in any predicted results. 

The CFD results for the 2-D engine configuration 
indicated the need for full 3-D analysis of the 
combustor flow field to properly predict the three- 
dimensional penetration and mixing inherent in 
the combustion process. An accurate modeling 
of the mixing process, especially for diffusion- 
controlled combustion, is required in order to 
adequately predict overall engine performance, 
with either a 1-D cycle code or a full-CFD analysis 
code. The forebody flow field boundary layer 
should also be included in any nose-to-tail 
analysis. The presence of this boundary layer 
can have a significant effect on the efficiency of 
the scramjet and on the performance of the entire 
vehicle. 

The analysis of the liquid oxygen-augmented 
preburning hydrocarbon scramjet indicated that 


the concept does indeed produce combustion 
of the main fuel within the scramjet engine. The 
preburning process provides a sufficiently- 
elevated temperature flow into the main fuel 
injector region to support immediate combustion 
of the gaseous ethylene fuel. However, 
because of the significant amount of unburned 
fuel at the combustor exit, there is a need for 
better mixing efficiency within the combustor. 
For the current fuel injector configuration, 
improved engine cycle performance could also 
be achieved at a lower overall engine 
equivalence ratio, limited by engine cooling 
requirements. 
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Fig. 1 Hypersonic vehicle concept with an 
integrated hydrocarbon scramjet 
engine. 





Fig. 2 Schematic diagram of free-shear layer. 



Fig. 3 Velocity profiles 
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Fig. 4 Turbulent kinetic energy profiles 



Fig. 5 Shear layer growth rate (vorticity 
thickness) 
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Fig. 9 Temperature contours 



Fig. 10 Temperature contours (Zoom view) 



Fig. 1 1 Fuel mass fraction contours 
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Fig. 15 Temperature contours at a) symmetry 
plane of top injectors, b) symmetry 
plane of bottom injectors. 



Fig. 13 Schematic diagram of the 3-D scramjet 
engine geometry 


Fig. 16 C 2 H 4 mass fraction contours at a) 
symmetry plane of top injectors, b) 
symmetry plane of bottom injectors. 



Fig. 14 Pressure contours at a) symmetry plane 
of top injectors, b) symmetry plane of 
bottom injectors. 
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Fig. 17 O 2 mass fraction contours at a) 
symmetry plane of top injectors, b) 
symmetry plane of bottom injectors. 
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Fig. 18 Velocity vectors in vicinity of preburner 
injector 



Fig. 19 Velocity vectors in vicinity of main 
burner injector 


x= 1.658m 



Figure 20: Temperature Contours at Indicated Axial Locations 
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Figure 21 : C2H4 Mass Fractions Contour at Indicated Axial Locations 


14 






Mass Flow Rate Concentrations 


x=1.658m 



Figure 22: H 2 O Mass Fraction Contours at Indicated Axial Locations 
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Fig. 23 C 2 H 4 and CO mass concentrations as a 
function of axial locatio 



X (m) 

Fig. 24 Temperature as a function of axial 
location 



Fig. 25 Pressure as a function of axial location 



Fig. 26 Axial velocity as a function of axial 
location 
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ABSTRACT 

The results of a feasibility study of a hypersonic 
wavehder research vehicle with a hydrocarbon 
scramjet engine are presented. The integrated 
waverider/scramjet geometry is first optimized 
with a vehicle synthesis code to produce a 
maximum product of the lift-to-drag ratio and the 
cycle specific impulse, hence cruise range. 
Computational fluid dynamics (CFO) is then 
employed to provide a nose-to-tail analysis of the 
system at the on-design conditions. Some 
differences are noted between the results of the 
two analysis techniques. A comparison of 
experimental, engineering analysis and CFO 
results on a waverider forebody are also included 
for validation. 

INTRODUCTION 

Interest in the development of various types of 
hypersonic vehicles has recently seen a 
resurgence. There exists an increasing need for 
hypersonic research vehicles (HRV) to 
demonstrate integrated aerodynamic, 
propulsion, and structural technologies for 
hypervelocity design and to develop a research 
database for reducing the risk involved in the 
development of operational hypersonic vehicles. 
Conceptual design activities are currently under 
way at NASA Ames Research Center to 
determine the feasibility of such a research 
vehicle utilizing near-term technology (Fig. 1). 
The objective of this research is to define an 
integrated hypersonic cruise vehicle that 
demonstrates sustained air-breathing hypersonic 
propulsion. 
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The research goals of this activity are to provide 
an understanding of the underlying physics, 
verification of design tools, and validation of the 
technologies and systems needed. The 
research requirements are classified into two 
main areas: 1) basic research, and 2) systems 
technology demonstration, which addresses 
programmatic research issues and overall vehicle 
system integration and performance. The 
disciplinary research requirements include 
aerodynamics and aero-thermodynamics of 
hypersonic flight, hypersonic air-breathing 
propulsion system performance, structures and 
materials characteristics, and finally 
instrumentation requirements for hypersonic 
vehicles. 

In order to accomplish these research goals, the 
mission of the HRV requires sustained cruise at 
various Mach numbers to address numerous 
hypersonic research requirements, including 
vehicle performance, real gas effects, boundary 
layer transition, shock-boundary layer interaction, 
turbulence modeling, propulsion integration, and 
structures / material performance. These 
requirements lead to the need for a high lift-to- 
drag ratio (L/D) to achieve the desired test times 
and to maximize the payload fraction to allow 
greater degrees of instrumentation. Waverider 
configurations have received a high degree of 
interest for their potentially high lift-to-drag ratios 
and their flow quality at the inlet plane [1]. These 
characteristics of waveriders are very desirable for 
cruise missions with integrated engines, hence 
HRV's. However, the (L/D) max of the 
propulsion-integrated configuration may be 
much lower than that of the pure waverider 
shape. 

For the HRV, most of the research requirements 
dictate a need for a hydrocarbon scramjet and/or 
ramjet operating between Mach numbers of 6 to 
8. For a Mach number range of 4 to 10, 
hydrocarbon fuels provide sufficient engine 
specific impulse (Isp) performance, heat sink 
capability, and offer the potential of reduced 
vehicle size compared with hydrogen-powered 
designs. In addition, the handling and 
infrastructure requirements for the hydrocarbon 
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fuels have a distinct advantage compared to 
cryogenic hydrogen. 

A common factor among all hypersonic vehicles 
is the need to effectively integrate the propulsion 
system with the airframe structure to maximize 
vehicle performance. The design must account 
for the aerodynamic heating, stability and control, 
and materials and structures. This poses a 
significant challenge to the designer of air- 
breathing aircraft since the problem becomes 
multidisciplinary. 

In this paper, the conceptual design and analysis 
of an air-breathing, Mach 8, waverider- 
configured, hypersonic testbed vehicle is 
performed. A comprehensive vehicle synthesis 
and design code is first used to define an optimal 
baseline configuration and CFD methods are 
then employed to refine the accuracy of the 
predicted performance. A nose-to-tail analysis is 
performed for the power-on flight condition at the 
design Mach number. 

CONCEPTUAL DESIGN TOOLS 

For the results presented later in this paper, two 
levels of analysis were performed. They include 
engineering analysis methods and computational 
fluid dynamics. The simplified engineering 
codes have been traditionally used in the 
conceptual design process to predict 
representative vehicle performance but have 
degraded accuracy in regions where the 
simplified assumptions break down. These 
regions can be numerous in a complex 
hypersonic vehicle design with integrated 
propulsion systems. Further, these simplified 
methods lack the capability to predict any 
unforeseen physics associated with a particular 
design. CFD, on the other hand, can significantly 
improve the accuracy and detail, but not without 
penalty. Significant computer resources can be 
required for a complete CFD analysis of the 
design. This section of the paper describes both 
analysis tools used in the overall design process. 

Hypersonic Vehicle Synt hesis Code 

The HAVOC hypersonic vehicle synthesis code 
[2] can be used to design, analyze and optimize 
hypersonic waverider configurations with or 
without an integrated hydrocarbon scramjet. The 
optimization methodology is that of Ref. 3. The 
aero/aerothermal and propulsion flow path 
techniques are briefly discussed below. 


The geometric definition of a hypersonic 
waverider configuration is computed by 
assuming that the lower surface of the waverider 
is a stream surface in an axisymmetric shock layer. 
Inverse design techniques are employed to 
determine this stream surface from a previously 
computed shock layer. The upper surface of 
the waverider is simply defined as the free-stream 
surface containing the waverider leading edge. 
The generating surface can be defined in either 
the free-stream (hence upper vehicle surface), 
or on the lower vehicle surface at an arbitrary 
longitudinal location. A sixth-order polynomial is 
used to describe the surface geometry. Solution 
of the real-gas Taylor-Macoll equations gives the 
inviscid flow properties throughout the shock 
layer. A simplified compressible boundary layer 
reference enthalpy method [4] is used to 
compute the local skin friction coefficient, which 
is used in turn to compute equilibrium radiation 
surface temperatures. No viscous-inviscid 
interactions are modeled in this engineering 
analysis approach. Leading-edge temperatures 
are computed using a swept-cylinder model [5]. 
Pressure lift and drag are computed by 
integration of the pressure coefficient over the 
surface of the vehicle. Base drag is computed 
using 70% vacuum pressure coefficient in the 
vehicle base region. 

The simplified, nose-to-tail propulsion model 
consists of an inviscid, 2-D, real gas, shock/weak- 
wave flow code coupled to a i-D 
subsonic/supersonic combustor analysis code. 
The shock/weak-wave code solves the inviscid 
inlet flow field as a function of vehicle 
forebody/ramp geometry, including cowl 
position, angle of attack and free stream Mach 
number. Equivalent 1-D flow properties are then 
computed at the inlet throat, and the 1-D 
combustor mass, momentum and energy 
equations with wall skin friction and heat transfer 
are solved stepwise through the burner. 
Combustor efficiency (i.e., heat release schedule 
as a function of combustor station) is input and 
was taken from Ref. 6 for the present results. 
The nozzle flow field is then computed from the 
combustor exit solution using the real-gas 
shock/weak-wave 2-D code, including nozzle 
and cowl flap geometry. First order estimates of 
axial and normal forces and pitching moment are 
thus computed as a function of vehicle geometry 
and flight condition. Overall propulsion system 
heat loads are then used to determine fuel inlet 
temperature or to compute required engine 
cooling equivalence ratio. 
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CFD Codes 

Computational fluid dynamics has and will 
continue to play an important role in the design 
and analysis of hypersonic systems. This is 
because ground based facilities are expensive to 
operate and in many cases cannot duplicate the 
exact flight conditions of such vehicles. 
Simplified methods such as that described above 
have limitations and cannot predict any 
unforeseen physics. Although studies have 
been performed that use CFD to investigate 
various components of hypersonic vehicles 
including waverider forebodies [7-9], complete 
power-on, nose-to-tail studies are scarce in the 
open literature. However with recent advances in 
both numerical algorithms and computer 
technology, such solutions are becoming 
possible. 

Requirements of a numerical algorithm are very 
demanding for hypersonic, nose-to-taii 
computations. The algorithm must be capable of 
predicting the three-dimensional flow of a highly 
turbulent mixture of reacting gasses with 
separated regions and strong flow field 
discontinuities. Multi-equation turbulence 
models become a requirement to accurately 
compute the numerous shear layers and 
boundary layers that are present. Further, 
upwind algorithms offer an appealing approach to 
solutions of hypersonic flow fields because of 
the ability to capture strong flow field 
discontinuities without user-specified smoothing 
terms. Finally, a strongly coupled, 
nonequilibrium chemistry model is required to 
compute the highly reactive combustion 
processes within a scramjet engine. 

For the present numerical analysis, the TUFF and 
STUFF codes of Ref. 10 were employed. The 
TUFF code is a time-marching code. It is 
generally used to obtain solutions in the 
subsonic or separated regions of the hypersonic 
flow field. Large run-times can prohibit the 
computation of an entire hypersonic flow field 
with a scheme of this nature. Hence, the second 
algorithm, STUFF, was developed. It employs a 
space-marching algorithm that can obtain a 
solution in relatively little computer time. For the 
nose-to-tail solution presented later, the STUFF 
code was employed to obtain nozzle and 
external solutions. The TUFF code was 
employed in Ref. 6 to obtain the solution of the 
current combustor. The combustor results in this 
paper were computed with the HAVOC code by 
employing the CFD-predicted heat release 
schedule of Ref. 6. 


The TUFF and STUFF codes offer many of the 
required features for accurate hypersonic flow 
field computations. Both codes employ 
nonequilibrium, equilibrium and perfect gas 
models. They employ a finite-volume philosophy 
to ensure that the schemes (including the 
boundary conditions) are fully conservative. 
Further, they obtain their upwind inviscid fluxes 
by employing a Riemann solver that fully 
accounts for the gas model used. This property 
allows the flow field discontinuities, such as 
shocks and shear layers, to be captured without 
significant amounts of smearing. Total variation 
diminishing (TVD) techniques are included to 
allow extension of the schemes to higher orders 
of accuracy without introducing* spurious 
oscillations. The schemes employ a strong 
coupling between the fluid dynamic and 
chemistry equations and are made fully implicit to 
eliminate the step-size restriction of explicit 
schemes. A fully conservative zonal scheme 
[11] has been implemented to allow solutions of 
geometrically complex problems. Turbulence 
models include both a zero [12] and two- 
equation [13] model with a correction for 
compressibility [14]. Finally, a sublayer 
approximation [15] is used in the space-marching 
algorithm to allow stable marching in the 
presence of a subsonic viscous layer. For the 
computations with hydrocarbon-air chemistry 
presented later, the thermodynamic data of Ref. 
16 and kinetics data of Ref. 17 were used. These 
codes have been validated on many classes of 
hypersonic flows including internal and external 
flow fields. 

For the waverider forebody validation case 
presented later in this paper, a modification to the 
viscosity model was required. Because of the 
very low temperatures experienced in the 
hypersonic test facilities, Sutherland's law for 
viscosity is no longer valid. The viscosity model 
of Ref. 18 was therefore employed in the CFD 
codes. This model was designed for hypersonic 
wind tunnel analysis and can be used to 
temperatures at which condensation begins. 
Also, base drag was computed assuming a 70% 
vacuum pressure coefficient, as in the 
engineering analysis above. 

Grid generation was accomplished by employing 
an interactive surface grid generator, S3D [19], 
and an interactive volume grid generator, 
HYPGEN-UI [20], for the external grids. Cross- 
section geometry obtained from the engineering 
analysis was directly used in the surface grid 
generation process. For the internal engine grid 
generation, an algebraic solver was employed. 
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The combination of these grid generation 
techniques proved to be quite effective in this 
current analysis. 

RESULTS 

Two sets of results are included in this paper: a 
pure waverider forebody test case and an 
integrated waverider/scramjet result. The first 
test case serves as a validation case for both the 
CFD and engineering analysis techniques on 
waverider geometries. The second case 
consists of the results of a conceptual design 
and analysis of a hypersonic waverider research 
vehicle with a hydrocarbon scramjet engine. 

Waverider Forebodv Result 

The first set of results presented here are for a 
pure waverider forebody at on-design 
conditions. The geometry and flow conditions 
are those defined by Reuss[21] at Ohio State 
University (OSU). In Ref. 21 a viscous-optimized 
waverider design was developed and tested. 
The MAX WARP code [22] was employed to 
generate a conical waverider geometry with a 
maximum lift-to-drag ratio at the OSU supersonic 
tunnel conditions. For this analysis, the free 
stream properties were those of the OSU 
supersonic facility and the geometry was 
restricted to a conical shock waverider resulting 
from a 12 degree cone. The freestream 
properties are listed in Table 1 . The optimization 
procedure accounted for both pressure and 
viscous forces including base drag. Throughout 
the optimization procedure, the coupling 
between the boundary layer and inviscid flow 
field was neglected. For this reason, the 
boundary layer displacement thickness was 
removed from the experimental model on both 
the upper and lower surfaces. The resulting 
waverider model was then installed in the OSU 
supersonic facility to provide surface pressure 
measurements only. The experimental model is 
0.1651 meters in length and is shown in Fig. 2. 


Property 

Value 

Mach 

8.00 

Reynolds No. 

1140000/m 

U» 

1227 m/s 

Poo 

.00377 kg/m 3 

Too 

58.6 K 

Gas 

Air 


Table 1. OSU tunnel operating conditions. 


The HAVOC code and both the TUFF and 
STUFF CFD codes were used to obtain flow field 
solutions of the OSU waverider at on-design 
conditions. The HAVOC results were obtained 
on the unaltered geometry since the effect of 
boundary layer displacement is neglected as in 
the MAXWARP code. However, the modified 
geometry was used in the CFD analysis since 
there is a direct coupling of the inviscid and 
viscous phenomena in a full numerical analysis. 

As mentioned earlier, the S3D and HYPGEN-UI 
codes were used for the grid generation 
procedure for the CFD results. Grid points were 
clustered near the leading edge of the waverider 
since high gradients can exist in that region. The 
rid spacing at the surface was set to 7.6 x 1 0* 
m. The outer boundary of the grid was placed 
well beyond the upper and lower bow shocks to a 
distance of 0.076m from the body. The surface 
grid and volume grid are shown in Fig. 3. The 
space-marching code, STUFF, was then used to 
obtain a solution about the entire geometry. This 
was possible since the nose was pointed 
allowing the Mach number in the marching 
direction to remain supersonic. This condition is 
required for stable marching. TUFF results were 
also obtained for comparison. Because the test 
time of the experiment was on the order of a 
minute, adiabatic boundary conditions were used 
in the CFD analyses. A perfect gas was assumed 
for these computations. 

Density contours of the CFD results at the last 
axial cross-section are shown in Fig. 4. This 
figure clearly shows the bow shock and boundary 
layer that are present in both the upper and 
lower flow fields. The boundary layers on both 
the top and bottom of the waverider are quite 
thick and comprise nearly one-third of the shock- 
layers. The effect of the relatively thick boundary 
layers on the inviscid flow field becomes 
apparent with the slightly detached shock wave 
from the waverider leading edge. Further, the 
inviscid portion of the flow field no longer exhibits 
a conical nature because of the influence of the 
boundary layer. 

The pressure contours at 95% of the body 
length are shown in Fig. 5 for all of the analysis 
techniques and the OSU experiment. All of the 
analysis results agree quite well at on-design 
conditions including TUFF and STUFF and 
engineering analysis. A grid refinement analysis 
was also performed to determine any grid 
dependency and very little was observed. The 
experimental results, on the other hand, show a 
14 percent higher pressure on the waverider 
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forebody than any of the analysis results. 
Comparison of STUFF results at a one-degree 
angle of attack agree much better with 
experiment. This seems to point to an 
uncertainty in the experimental angle of attack. 
For the experimental results presented here, a 
one-degree error in the experimental angle of 
attack was possible and fails within the 
experimental uncertainty range [21]. Further 
experiments are planned to resolve this 
discrepancy. 

Figure 5 also shows the enhanced ability of the 
CFD techniques to predict viscous, waverider 
forebody flow fields. This is apparent in the large 
pressure rise predicted at the waverider leading 
edge by the CFD techniques. This pressure rise 
is caused by compression waves emanating from 
the rapidly growing boundary layer in the vicinity 
of the leading edge. This effect is entirely 
neglected in the engineering code results since 
no coupling is allowed between the viscous and 
inviscid analysis. This effect, however, is also 
present in the experimental results. 



MAXWARP 

HAVOC 

STUFF 

Cl 

0.0870 

0.0885 

0.0844 

C D 

0.0241 

0.0257 

0.0313 

L/D 

3.62 

3.45 

2.69 


Table 2. Comparison of predicted aerodynamic 
coefficients. 

Table 2 shows the predicted aerodynamic 
coefficients of both engineering codes and of 
the CFD analysis. Since only very little difference 
was observed between the TUFF and STUFF 
aerodynamic coefficients, only the STUFF 
results are included in this table. All of the 
predicted lift coefficients agree quite well even 
though the leading edge pressure rise was 
present in the CFD results. The increased lift 
caused by the pressure rise on the lower surface 
is counteracted by a similar pressure rise on the 
upper surface in the CFD results. The drag 
coefficient predicted by the CFD analysis 
however is slightly larger. This is simply 
explained by increased pressure drag on the 
forebody resulting from the boundary layer 
displacement especially at the leading edge. 


Hypersonic Waverider Research Vehicle Design 

and Analysis 

Design Optimization 

As opposed to accelerator-type missions (e.g., 
SSTO) where the mass capture characteristics of 
the vehicle are most important, cruise 
configurations place an emphasis on the 
aerodynamic performance of the design. Hence 
vehicle lift-to-drag ratio (L/D) or the product of 
L/D and engine Isp (a parameter proportional to 
Brequet range factor) becomes more important. 
A waverider configuration, with high hypersonic 
L/D potential, was selected as the baseline 
configuration for the present study of the HRV. 
The configuration/engine installation was 
numerically optimized to maximize (L/D) x Isp, 
using forebody shape, ramp angles and cowl 
position as optimization parameters. For this 
analysis, the forebody is a Mach 8 waverider 
configuration and the ramp angles and cowl 
position were designed to produce first and 
second ramp shock-on-lip and cowl shock-on- 
shouider. The engine geometry and operating 
parameters were held fixed in the optimization 
procedure. The numerical optimization was 
performed using the HAVOC cycle code. 

The design parameters and constraint functions 
used for optimization are listed below. There 
were 15 design parameters used in the 
optimization process. Six parameters defined 
the waverider generating surface (hence the 
vehicle shape); eight parameters defined the 
ramps, cowl, nozzle geometries, and leading- 
edge radius; the last design parameter was free 
stream dynamic pressure. Eight constraints were 
used in the optimization process, including 
engine throttle greater than or equal to the 
required cooling equivalence ratio, leading edge 
equilibrium radiation temperature limited to 
3200°F, wing tip closure angle limited to 10 
degrees, and the vehicle structural thickness limit 
at vehicle aft end, accounting for nozzle 
integration. The volume was constrained to 
1.7% of vehicle length cubed. Finally, the 
vehicle width to length ratio was less than or 
equal to 0.75. 

For the Mach 8 design, a generating shock angle 
of 12 degrees was arbitrarily selected. The 
waverider shape optimization process involved 
planform shape changes to sweep the leading 
edge in order to alleviate high heating rates at the 
higher free stream dynamic pressures, traded off 
against leading edge radius and associated 
leading edge bluntness drag. The design 
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optimization process of the inlet (ramp positions 
and ramp angles) resulted in the two ramp shocks 
converging on the cowl lip (shock-on-lip) then 
reflected to the shoulder of the combustor 
entrance (shock-on-shoulder). The geometric 
contraction ratio was approximately 14, with a 
resulting pressure at the combustor inlet of 
about one atmosphere. Preliminary 
performance estimates indicated a required 
vehicle length of 23 ft, with an overall vehicle 
body density of 20 lbs/ft**3. A Jift-to-drag ratio of 
4.3 was achieved, with a cowl-to-taii Isp of 746 
seconds at an assumed combustor efficiency of 
95% and at an equivalence ratio of 1 .0. with a 
free stream dynamic pressure of 900 psf, 
resulting in an initial cruise altitude of 92,500 ft. 
An engine width of roughly 0.762m. produced a 
net thrust that is equal to the net drag of the 
vehicle at the design point. For the results 
presented below, the engine combustion 
efficiency of 55% was taken from Ref. 6 

Scramjet Engine Concept 

The slow reaction rates of a hydrocarbon/air 
mixture in a supersonic stream can have a 
significant impact on the inlet/combustor design. 
Because of the size limitation of an engine on a 
small research vehicle, a mechanism is required 
to provide sufficient fuel/air temperatures for 
burning within the combustor. The design 
concept of Ref. 6 is employed here and uses a 
liquid oxygen (LOX) augmented pre-burner 
located upstream of the main fuel injectors to 
promote burning in the combustor. Because the 
LOX would be stored onboard the research 
vehicle, the required pre-burning should be 
kept to a minimum to reduce the impact on the 
vehicle design and gross weight. For the 
relatively-short hypersonic research mission 
cruise times (5 to 10 minutes), this additional on- 
board mass should not have a significant impact 
on the HRV design. 

Preliminary HAVOC design results at a Mach 
number of 8 for the waverider HRV indicated that 
the air throat temperature was 833°K. Given that 
the spontaneous ignition temperature of an 
ethylene/oxygen mixture is roughly 780°K[23], a 
guideline of 1000°K was used as the design 
temperature at the mainbumer station to insure 
that combustion was indeed present. The 
required prebumer fuel flow and prebumer O/F 
ratio were then computed. To achieve the 
design temperature at the main fuel injector 
station, approximately 2.5% of the overall 
stoichiometric engine fuel flow with a prebumer 
O/F ratio of 1 .71 (twice prebumer stoichiometric) 


is required. The CFD results of Ref. 6 confirmed 
that burning was present at the main burner 
station under these conditions. A heat balance 
on the vehicle and engine was used to compute 
the fuel total temperature at injection. Detailed 
engine operating conditions and geometry are 
presented in Ref. 6. 

Nose-to-Tail Analysis 

Nose-to-tail analyses utilizing both CFD and cycle 
codes were performed. The geometries for both 
analyses, however, were somewhat different. 
The ramp sidewalls in the cycle analysis were 
assumed to have zero thickness and were 
aligned with the forebody flow. They therefore 
resulted in zero pressure drag with only a nominal 
addition to the viscous drag. The CFD geometry, 
on the other hand, contained ramp sidewalls with 
an included angle of 21 deg. Because the 
sidewall attaches to the cowl lip and because the 
cowl protruded far below the lower surface of the 
waverider, the sidewall contained a significant 
amount of additional volume. This added volume 
was exposed to the forebody flow and produced 
additional drag to the vehicle. 

The grid generation for the CFD nose-to-tail 
analysis is described below. The grid generation 
tools used for the external grids are the S3D and 
HYPGEN-UI codes mentioned previously. The 
external portion was divided into three axial 
sections: nose-to-inlet, inlet-to-exit, and exit-to- 
end sections. Grid generation was performed 
separately on each of these portions. For the 
surface grid, points were clustered on the 
waverider leading edge, on the sidewall leading 
edge and on all of the other convex comers of 
the vehicle. The lower surface grid is shown in 
Fig. 6. The outer boundary of the volume grid 
was placed well beyond the anticipated bow 
shock to a distance of roughly three meters from 
the vehicle surface. The spacing of the first point 
from the surface was set to 1 .0 x 1 0' 5 m. The grid 
dimensions for the three external grids measured 
102x70x144, 60x70x129 and 12x70x130 from 
nose to tail. The values of these dimensions 
correspond to the number of grid points in the 
streamwise direction, circumferential direction, 
and radial direction. 

For the internal grids, the engine was divided into 
three separate flow paths starting at the inlet 
face. These flow paths then merged at the 
internal nozzle entrance. An algebraic grid 
generation routine was used to generate ail of 
the interior grids. The internal grids measured 
7x79x50 for the inboard engine, 7x79x99 for the 
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outboard engines and 3x79x79 for the internal 
nozzle section. The number of grid points in the 
streamwise direction is small since the space- 
marching scheme interpolates for additional grid 
planes as they are needed. 

The space-marching scheme, STUFF, was 
employed to obtain the nozzle and external CFD 
results for the waverider HRV. A space-marching 
solution is possible if the flow field is void of 
streamwise subsonic and separated areas. The 
space-marching solution began at the nose of 
the vehicle by setting the dependent variables in 
all of the cells to be freestream. The solution 
progressed by marching downstream through 
each of the grids until the aft end of the vehicle 
was reached. A fully conservative patching 
scheme [11] was employed to transfer the 
solution from one grid to the next. The frozen 
chemistry option was employed for the CFD 
results including the nozzle portion. This was 
sufficient since combustion is slowed 
considerably at the combustor exit [6]. 

For the nose-to-tail results presented here, the 
HAVOC cycle code was used to provide the 
solution within the inlets and combustor. This 
was done to avoid the computational expense of 
obtaining time-marching solutions of both the 
outboard and inboard engines with 48 fuel 
injection ports in each. A time-marching solution 
is required in portions of the engine because the 
presence of axially-subsonic and axially- 
separated regions prohibits space-marching. 
The CFD forebody results were averaged at the 
inlet face to provide one-dimensional inlet 
conditions for the cycle code. The presence of 
an oblique cowl shock was accounted for and 
provided the mechanism to turn the flow parallel 
to the engine cowl. The CFD-predicted heat 
release schedule of Ref. 6 was used in the one- 
dimensional cycle analysis of the combustor. 
CFD then picked up again at the exit of the 
combustor. 

The two-dimensional, shock/weak wave results 
of the HAVOC analysis are plotted in Fig. 7. This 
figure shows the shock, contact surface, and 
expansion wave arrangement for the current 
integrated design on the vehicle symmetry 
plane. The intent of the forebody/ ramp design is 
ciearty shown in this figure. Both of the ramp 
shocks impinge on the cowl leading edge and 
the cowl shock impinges on the shoulder. The 
bow shock, however, lies outside of the inlet. 
This is because at Mach 8, with a fixed throat 
height, having the bow shock on the cowl lip can 
produce shock-on-ramp with possible boundary 


layer separation and engine unstart. The nozzle 
geometry was restricted to two planar sections 
and the corresponding expansion and shock are 
shown in Fig. 7. 

Figure 8 is a plot of the surface pressure on the 
keel line of the vehicle. Both the CFD and cycle 
results are plotted. Agreement is excellent 
except in the vicinity of inviscid/viscous 
interactions. Comparison of the pressure on the 
waverider forebody are within 2 percent except at 
the leading edge where CFD predicts a pressure 
spike resulting from boundary layer 
displacement. The ramp pressures also agree 
quite well except for the asymptotic behavior of 
the CFD results that is typical of shock/boundary 
layer interactions. Other reasons for this slight 
discrepancy are three-dimensional effects. One 
such effect is encountered as the planar ramp 
emerges from the curved waverider forebody. 
The intersection line bends downstream from the 
vehicle symmetry plane. This three-dimensional 
effect reduces the pressure near the intersection 
region. Figure 8 also shows very good 
agreement on the outer surface of the cowl and 
on the nozzle surface. Comparison of surface 
pressure on the upper surface showed that the 
CFD predicted pressures were 3-10% higher 
than the engineering predictions. This is also 
explained by viscous/inviscid interactions. 

The pressure contours in Fig. 9 show the shock 
and expansion waves that are predicted with 
CFD. The location and strength of these waves 
are in very good agreement with those predicted 
by the cycle code (Fig. 7). The CFD results, 
however, predict that the ramp shocks lie slightly 
outside of the inlet instead of on the cowl lip. 
This difference is attributed to shock/boundary 
layer interactions and to three-dimensional 
effects. The boundary layers on the upper and 
lower surface are visible in the density contours 
of Fig. 10. 

Figures 11 and 12 show the CFD predicted 
pressure and density contours, respectively, on 
crossflow planes at various axial stations. These 
stations correspond to the following locations: a) 
beginning of the second ramp, b) at the inlet, c) 
at the middle of the cowl, d) at the engine exit 
and e) at the vehicle end. The two ramp shocks 
are visible in the pressure contours of Fig. 11a 
and clearly exhibit a 3-dimensional behavior. The 
bow shock is attached to the leading edge of the 
waverider and produces a very clean flow field on 
the forebody absent of the ramp geometry. 
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Figures 11b and 12b show that the flow into the 
inlet is fairly clean except for a weak shock that 
emanates from the ramp sidewall and interacts 
with the forebody boundary layer. Also evident 
in Fig. 12b is a thickening of the ramp boundary 
layer thickness nearer the ramp sidewall. This is 
caused by the tendency of the flow to spill off of 
the first ramp resulting in a sideways velocity 
component. This effect washes the boundary 
layer away from the symmetry plane of the 
vehicle. It is then stopped before the inlet station 
by the sidewalls on the second ramp. Another 
feature that should be noted in Figs. 11b and 
12b is the shock emanating from the outer 
surface of the ramp sidewall. This feature is 
absent in the cycle analysis since the sidewalls 
were assumed to have zero thickness. The 
resulting high pressure region adds considerably 
to the net drag of the vehicle since the outer 
sidewall has an area component in the 
streamwise direction. This effect resulted in a 
negative net-thrust predicted by the nose-to-tail 
analysis. This therefore leads to a need to 
reduce the thickness of the sidewalls. The 
pressure and density contours on the external 
portion of the engine at the midpoint of the cowl 
are plotted in Figs 11c and 12c. The ramp 
shocks have clearly merged with the bow shock 
at this point and an expansion fan emanating 
from the lower surface of the cowl is evident. 

Figures 11(d-e) and I2(d-e) depict the CFD 
predicted flow field on the aft portion of the 
vehicle. This flow field region is very complex 
due to the presence of the engine, although a 
significant portion of the waverider flow field 
remains unaffected. The bow shock remained 
attached to the leading edge for the entire length 
of the vehicle. This feature is desirable since any 
spillage of the high pressure gasses onto the 
upper surface would reduce the performance of 
the vehicle. Two flow field features that are 
cfearty visible in the nozzle portion of the flow are 
the initial nozzle expansion and the shock 
caused by the second nozzle plane. These 
features appear to be nearly two-dimensional, 
leading to the good agreement with the cycle 
analysis (Fig. 8). 

Figure 13 shows the pressure contours on the 
lower surface of the vehicle. The impact on the 
surface pressure by the presence of the 
integrated scramjet is visible in this figure. The 
lower surface is comprised of a large region of 
undisturbed flow with a fairly constant pressure. 
The remaining surface is exposed to the 
numerous shocks and expansion fans that 
originate from the integrated propulsion system. 


These include the ramp shocks, the sidewall 
shocks, the exterior cowl expansions, the engine 
shroud expansion, and the numerous nozzle 
waves. 

Of particular interest in the current design is the 
forebody ramp system and the resulting flow 
field. Due to the curvature of the waverider 
forebody, an oblique intersection results 
between it and the first ramp. This intersection is 
visible in the pressure contours of Fig. 13. The 
strong pressure gradients that are present in this 
vicinity cause the flow to be diverted away from 
the vehicle symmetry plane beginning at the 
origin of the first ramp. This effect tends to 
reduce the amount of mass provided to the inlet. 
An intersection that is less oblique would reduce 
the spillage since the pressure gradient would be 
aligned more with the flow direction. This can be 
accomplished by reducing the curvature of the 
waverider forebody at the symmetry plane. In the 
current design, this spillage was effectively 
halted at the second ramp by using sidewalls. 
This is shown in the surface streamlines of Fig. 
14. 

CONCLUSIONS 

Two methods were demonstrated for the 
conceptual design and analysis of a hypersonic 
waverider research vehicle: an engineering 
analysis code with a simplified nose-to-tail flow 
field analysis capability and a complete 3-D CFD 
code with a hydrocarbon/air capability. The 
methods have been shown to agree well for 
waverider forebodies except in the vicinity of the 
leading edge where CFD exhibits a superior 
capability to predict the strong viscous/inviscid 
interaction. Comparison with the OSU waverider 
experimental results were inconclusive since the 
angle of attack of the experiment was 
questionable. The coupling of the cycle code 
with CFD to compute a single flow field was also 
demonstrated. Even though the cycle code has 
a significantly degraded ability to predict detail, 
this coupling proved to be useful because of the 
reduced computer resources required. 

The nose-to-tail analysis of the waverider HRV 
clearly showed benefits of the current design 
and also revealed areas for improvement. The 
waverider forebody combined with the current 
ramp system provided a uniform flow to the inlet 
of the scramjet. Further, the forebody shock 
remained attached to the leading edge for the 
entire length of the vehicle. This avoided 
spillage of the high pressure air onto the upper 
surface which could significantly reduce lift. The 
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analysis also showed that the ramp sidewalls 
need to be thinner since they protrude far into 
the forebody flow field and cause unnecessary 
drag to the vehicle. 

FUTURE RESEARCH 

Future research includes addressing the issues 
found in the analysis of the HRV design. These 
include the engine performance, nozzle 
performance, side-wall drag and inlet efficiency. 
Further studies also include structural 
analysis/weight estimation, power-on/power-off 
stability and control analysis and off-design 
performance studies. 

ACKNOWLEDGMENTS 

This work was partially supported by NASA Ames 
Research Center under grants NCC2-498 and 
NCC2-746. The authors would like to 
acknowledge Dr. Hirokazu Miura at NASA Ames 
Research Center for providing a waverider model 
deformation analysis and Mr. Thomas Whitaker of 
Sterling Software for graphics and technical 
support. 

REFERENCES 

[1 ] O'Neill, M. K., and Lewis, M. J., "Optimized 
Scramjet Integration on a Waverider," AIAA 
Paper 91-1693, June 1991. 

[2] Bowles. J. V. Conceptual Studies 
Activities. Proceedings of the Second 
National Aerospace Plane Symposium. 
Applied Physics Laboratory, Laurel, MD, 
Nov. 1986. 

[3] Vanderplaats, G. N. and Hansen, S. 
R.,"DOC Users Manual", Version 1.00, 
VMA Engineering, Goleta, CA, 1989. 

[4] Schlichting, H., "Boundary-Layer Theory", 
Seventh Edition, McGraw-Hill Book Co., 
New York, 1979. 

[5] Anderson, J. D., "Hypersonic and High 
Temperature Gas Dynamics," First Edition, 
McGraw-Hill Book Co., New York, 1989. 

[6] Molvik, G. A., Bowles, J.V. and Huynh, L. 
C., "Analysis of a Hydrocarbon Scramjet 
with Augmented Presuming," AIAA Paper 
92-3425, July, 1992. 

[7] Jones, K. D., "Application of a Supersonic 
Full Potential Method for Analysis of 


Waverider Configurations,” NASA TP- 
2608, 1986. 

[8] Takashima, N., and Lewis, M. J., "Navier 
Stokes Computation of a Viscous 
Optimized Waverider,” AIAA Paper 92- 
0305, Jan. 1992. 

[9] Jones, K. D., Bauer, S. X. S., and 
Dougherty, F. C., "Hypersonic Waverider 
Analysis: A Comparison of Numerical and 
Experimental Results,” AIAA Paper 91- 
1696, June 1991. 

[10] Molvik, G. A. and Merkle, C. L., "A Set of 
Strongly- Coupled Upwind Algorithms for 
Computing Flows in Chemical 
Nonequilibrium," AIAA Paper 89-0199. 
Jan. 1989. 

[11] Klopfer, G. H. and Molvik G. A., 
"Conservative Multizonal Interface 
Algorithm for the 3-D Navier-Stokes 
Equations," AIAA Paper 91-1601, June 
1991. 

[12] Baldwin, B. S. and Lomax, H., “Thin Layer 
Approximation and Algebraic Model for 
Separated Turbulent Flows.” AIAA Paper 
78-257, Jan. 1978. 

[13] Jones, W. P. and Launder, B. E., "The 
Prediction of Laminarization with a Two- 
Equation Model of Turbulence," 
International Journal of Heat and Mass 
Transfer, Vol. 15, 1972, pp. 303-314. 

[14] Zeman, O., "Compressible Turbulence 
Subjected to Shear and Rapid 
Compression," Eighth Symposium on 
Turbulent Shear Flows, Munich, Germany, 
Sept, 1991. 

[15] Vigneron, Y. C., Rakich, J. C. and 
Tannehill, J. C., "Calculation of supersonic 
Viscous Flow over Delta Wings with Sharp 
Subsonic Leading Edges," AIAA Paper 
78-1137, July 1978. 

[16] Gordon, S. and McBride, B. J., "Computer 
Program for Calculation of Complex 
Chemical Equilibrium Compositions, 
Rocket Performance, Incident and 
Reflected Shocks, and Chapman-Jouguet 
Detonations," NASA SP-273, 1976. 

[17] Westbrook, C. K. and Dryer, F. L., 
"Simplified Reaction Mechanisms for the 


9 



Oxidation of Hydrocarbon Fuels in 
Flames,", Combustion Science and 
Technology, Vol. 27, 1981, pp. 31-43. 

[18] Reinecke, W. G., "Charts for Use with 
Hypersonic Air Wind Tunnels," Aerospace 
Research Laboratories Report, ARL 64-56, 
April, 1964. 

[19] Luh, R. C., Pierce, L. E., Yip, D., 
"Interactive Surface Grid Generation," AIAA 
Paper 91-0796, Jan. 1991. 

[20] Chan, W. M. and Steger, J. L., 
"Enhancements of a Three-Dimensional 
Hyperbolic Grid Generation Scheme," 
Applied Mathematics and Computation, 
Vol 51, pp. 181-205, 1992. 

[21] Reuss, R. L., "Experimental Study of 
Hypersonic Waveriders," Masters Thesis, 
Ohio State University, 1993. 

[22] Corda, S., and Anderson, J. D. Jr., 
"Viscous Optimized Hypersonic 
Waveriders Designed from Axisymmetric 
Flow Fields." AIAA Paper 88-0369, Jan. 
1988. 

[23] Barnett, H. C. and Hibbard, R. R., "Basic 
Considerations in the Combustion of 
Hydrocarbon Fuels with Air," NACA Report 
1300, 1959. 



FIGURES 



Fig. 1 Hypersonic vehicle concept with an 
integrated hydrocarbon scramjet 
engine. 



Fig. 3 Grid used for the CFD analysis of the 
OSU waverider. 
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Fig. 4 CFD predicted density contours at the 
aft end of the OSU waverider. 



Fig. 5 Pressure comparison at x=0.1 56845m 



Fig. 7 HAVOC predicted flow field on the 
lower symmetry plane. 




Fig. 9 CFD predicted pressure contours on 
vehicle symmetry plane. 




Fig. 6 Lower surface grid for the current 
waverider HRV design. 


Fig. 10 CFD predicted density contours on 
vehicle symmetry plane. 


1 1 




Fig. 1 1 CFD predicted pressure contours at 
various axial locations. 


Fig. 12 CFD predicted density contours at 
various axial locations. 







Fig. 13 CFD predicted surface pressure plot. 



Fig. 14 CFD predicted surface streamlines on 
the forebody ramp system. 
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ABSTRACT 

The results of a feasibility study of a hypersonic 
waverider research vehicle with a hydrocarbon 
scramjet engine are presented. The scramjet 
engine concept consists of a preburner into 
which a small amount of fuel is burned with on- 
board liquid oxygen and injected into the airflow, 
upstream of the main fuel injector locations, thus 
ensuring that main fuel combustion is present 
and. uninterrupted. The integrated 
wave rider/sc ramjet geometry is optimized with a 
vehicle synthesis code to produce a maximum 
product of the lift-to-drag ratio and the cycle 
specific impulse, hence cruise range. 
Computational fluid dynamics is employed to 
provide engine performance and a nose-to-tail 
analysis of the vehicle at the on-design 
conditions. Comparisons are made between the 
results of the two analysis techniques and some 
differences are noted. 

INTRODUCTION 

Interest in the development of various types of 
hypersonic vehicles has recently seen a 
resurgence. There exists an increasing need for 
hypersonic research vehicles (HRV) to 
demonstrate integrated aerodynamic, 
propulsion, and structural technologies for 
hypervelocity design and to develop a research 
database for reducing the risk involved in the 
development of operational hypersonic vehicles. 
Conceptual design activities are currently under 
way at NASA Ames Research Center to 
determine the feasibility of such a research 
vehicle utilizing near-term technology (Fig. 1). 
The objective of this research is to define an 
integrated hypersonic cruise vehicle that 
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demonstrates sustained air-breathing hypersonic 
propulsion. 

The research goals of this activity are to provide 
an understanding of the underlying physics, 
verification of design tools, and validation of the 
technologies and systems needed. The 
research requirements are classified into two 
main areas: 1) basic research, and 2) systems 
technology demonstration, which addresses 
programmatic research issues and overall vehicle 
system integration and performance. The 
disciplinary research requirements include 
aerodynamics and aero-thermodynamics of 
hypersonic flight, hypersonic air-breathing 
propulsion system performance, structures and 
materials characteristics, and finally 
instrumentation requirements for hypersonic 
vehicles. 

In order to accomplish these research goals, the 
mission of the HRV requires sustained cruise at 
various Mach numbers to address numerous 
hypersonic research requirements, including 
vehicle performance, real gas effects, boundary 
layer transition, shock-boundary layer interaction, 
turbulence modeling, propulsion integration, and 
structures/material performance. These 
requirements lead to the need for a high lift-to- 
drag ratio (UD) to achieve the desired test times 
and to maximize the payload fraction to allow 
greater degrees of instrumentation. Waverider 
configurations have received a high degree of 
interest for their potentially high lift-to-drag ratios 
and their flow quality at the inlet plane [1 j. These 
characteristics of waveriders are very desirable for 
cruise missions with integrated engines, hence 
HRVs. However, the (L/D) ma x the propulsion- 
integrated configuration may be much lower than 
that of the pure waverider shape. 

For the HRV, most of the research requirements 
dictate a need for a hydrocarbon scramjet and/or 
ramjet operating between Mach numbers of 6 to 
8. For a Mach number range of 4 to 10, 
hydrocarbon fuels provide sufficient engine 
specific impulse (Isp) performance, heat sink 
capability, and offer the potential of reduced 
vehicle size compared with hydrogen-powered 
designs. In addition, the handling and 
infrastructure requirements for the hydrocarbon 
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fuels have a distinct advantage compared to 
cryogenic hydrogen. 

The slow reaction rates of a hydrocarbon/air 
mixture in a supersonic stream can have a 
significant impact on the iniet/combustor design. 
Because of the size limitation of an engine on a 
small research vehicle, a mechanism is required 
to provide sufficient fuel/air temperatures for 
burning within the combustor. I he concept of 
employing an in-stream, embedded ramjet as a 
pilot light has been proposed [2-3] and appears 
to be a promising technique for maintaining 
combustion. An alternative is to use a liquid 
oxygen (LOX) augmented pre-burner located 
upstream of the main fuel injectors to promote 
burning in the combustor. Because the LOX 
would be stored onboard the research vehicle, 
the required pre-burning should be kept to a 
minimum to reduce the impact on the vehicle 
design and gross weight. For the relatively short 
hypersonic research mission cruise times (5 to 10 
minutes), the impact of this additional onboard 
mass will not have a significant impact on the 
HRV design. 

A common factor among all hypersonic vehicles 
is the need to effectively integrate the propulsion 
system with the airframe structure to maximize 
vehicle performance. The design must account 
for the aerodynamic heating, stability and control, 
and materials and structures. This poses a 
significant challenge to the designer of air- 
breathing aircraft since the problem becomes 
multidisciplinary. 

In this paper, the conceptual design and analysis 
of an air-breathing, Mach 8, waverider- 
configured, hypersonic testbed vehicle is 
performed. The feasibility of a hydrocarbon 
scramjet engine utilizing a LOX/preburner 
concept is first addressed[4] and then a 
comprehensive vehicle synthesis and design 
code is used to define an optimal baseline 
configuration^]. Computational Fluid Dynamic 
(CFD) methods were employed throughout the 
design process to refine the accuracy of the 
predicted performance. CFD results are included 
of the current scramjet engine concept and of the 
integrated vehicle at on design conditions. 

CONCEPTUAL DESIGN TOOLS 

For the results presented later in this paper, two 
levels of analysis were performed including: 
engineering analysis and computational fluid 
dynamics. Engineering codes have been 
traditionally used in the conceptual design 


process to predict representative vehicle 
performance but have degraded accuracy in 
regions where simplifying assumptions break 
down. These regions can be numerous in a 
complex hypersonic vehicle design with 
integrated propulsion systems. Further, these 
simplified methods lack the capability to predict 
any unforeseen physics associated with a 
particular design. CFD, on the other hand, can 
significantly improve the accuracy and detail, but 
not without penalty. Significant computer 
resources can be required for a complete CFD 
analysis of the design. This section of the paper 
describes both analysis tools used in the present 
design process. 

Hypersonic Vehicle Synthesis Coda 

The HAVOC hypersonic vehicle synthesis code 
[6] can be used to design, analyze and optimize a 
hypersonic waverider configuration including an 
integrated scramjet engine. The optimization 
methodology utilized in the HAVOC code is 
detailed in Ref. 7. The aero/aerothermai and 
propulsion flow path techniques are briefly 
discussed below. 

The geometric definition of a hypersonic 
waverider configuration is computed by 
assuming that the lower surface of the waverider 
is a stream surface in an axisymmetric shock layer. 
Inverse design techniques are employed to 
determine this stream surface from a previously 
computed shock layer. The upper surface of the 
waverider is simply defined as the free-stream 
surface containing the waverider leading edge. 
The generating surface can be defined in either 
the free-stream (hence upper vehicle surface), or 
on the lower vehicle surface at an arbitrary 
longitudinal location. A sixth-order polynomial is 
used to describe the surface geometry. Solution 
of the real-gas Taylor-Macoil equations give the 
inviscid flow properties throughout the shock 
layer. A simplified compressible boundary layer 
reference enthalpy method [8] is used to 
compute the local skin friction coefficient, which 
is used in turn to compute equilibrium radiation 
surface temperatures. No viscous-inviscid 
interactions are modeled in this engineering 
analysis approach. Leading-edge temperatures 
are computed using a swept-cylinder model [9]. 
Pressure lift and drag are computed by 
integration of the pressure coefficient over the 
surface of the vehicle. Base drag is computed 
using 70% vacuum pressure coefficient in the 
vehicle base region. 
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The simplified, nose-to-tail propulsion model 
consists of an inviscid, 2-D, real gas, shock/weak- 
wave flow code coupled to a 1-D 
subsonlc/supersonic combustor analysis code. 
The snock/weak-wave code solves the inviscid 
inlet flow field as a function of vehicle 
forebody/ramp geometry, including cowl 
position, angle of attack and free stream Mach 
number. Equivalent 1-D flow properties are then 
computed at the inlet throat, and the 1-D 
combustor mass, momentum and energy 
equations with wall skin friction and heat transfer 
are solved stepwise through the burner. 
Combustor efficiency (i.e., heat release schedule 
as a function of combustor station) was taken 
from the engine CFD results for the present 
analysis. The nozzle flow field is then computed 
from-the combustor exit solution using the real- 
gas shock/weak-wave 2-D code, including nozzle 
and cowl flap geometry. First order estimates of 
axial and normal forces and pitching moment are 
thus computed as a function of vehicle geometry 
and flight condition. Overall propulsion system 
heat loads are then used to determine fuel inlet 
temperature or to compute required engine 
cooling equivalence ratio. 

CFD Codes 

For the present numerical analysis, the TUFF and 
STUFF codes of Ref. 10 were used since they 
offer many of the features required for accurate 
hypersonic flow field computations including 
upwind fluxes and fully coupled chemistry. The 
TUFF code is a time-marching code. It is 
generally used to obtain solutions in the 
subsonic or separated regions of the hypersonic 
flow field. The STUFF code employs a space- 
marching algorithm that can obtain a solution in 
relatively little computer time. For the results 
presented later, the STUFF code was employed 
to obtain nozzle and external solutions. The 
TUFF code was used to obtain the solution within 
the combustor. 

Simplified reaction mechanisms for the 
hydrocarbon combustion process offer an 
appealing alternative to exhaustive computations 
with large chemistry systems. In this approach, 
various species and reactions are combined and 
simplified while preserving the net effect of the 
reaction processes. Since kinetic details are not 
required in the present study, a simplified 
reaction mechanism was sufficient. For the 
scramjet propulsion system analysis, the two-step 
reaction mechanism of Westbrook and Dryer [1 1 ] 
was employed to address the combustion of fuel 
with air. For hydrocarbon scramjet propulsion, 


the liquid fuel can be used as a coolant on various 
portions of the aircraft. This results in an 
endothermic reaction that can vaporize and 
dissociate the liquid fuel. Gaseous ethylene was 
used in this analysis, as a surrogate fuel intended 
to represent the products of this endothermic 
reaction. The thermodynamic and transport 
properties for the individual species in the 
hydrocarbon/air mixture were obtained from Ref. 
12 . 

Modeling scramjet flow fields with CFD requires 
an advanced turbulence model capable of 
accurately accounting for compressible turbulent 
shear layers and jets. This was accomplished in 
the present study with the incorporation of the 
low Reynolds number K-e turbulence model 
originally developed by Jones and Launder [13j. 
The compressibility correction of Zeman (1 4] was 
also included in the two-equation formulation to 
improve the computation of compressible shear 
layers. The turbulence model was transformed 
to a generalized finite-volume coordinate system 
and strongly coupled with the existing flow 
solver, including the viscous and inviscid flux 
computation and the source term treatment. 

Grid generation was accomplished by employing 
an interactive surface grid generator, S3D [15], 
and an interactive volume grid generator, 
HYPGEN-UI [16], for the external grids. Cross- 
section details determined by the HAVOC code 
were used in the surface grid generation 
process. For the internal engine grid generation, 
an algebraic solver was employed. The 
combination of these grid generation techniques 
proved to be quite effective and timely. 


RESULTS 

Two sets of results are included in this paper. ' 
The first set of results addresses the design and 
performance of a single hydrocarbon scramjet 
utilizing an augmented preburner upstream of 
the main fuel injectors. The second set consists 
of the results of a conceptual design and analysis 
of a hypersonic waverider research 
vehicle(HWRV) with a hydrocarbon scramjet 
engine. The conceptual design of the HWRV 
relied heavily on the CFD results of the scramjet 
engine. 

Hydrocarb on Scramiet Result 

The initial geometric definition for the scramjet 
engine, including throat height, shock isolator 
length, and combustor length, and combustor 
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area ratio was taken from Ref. 3. The embedded 
ramjet section was removed and a 
hydrocarbon/LOX preburner was added. A 
mixing section aft of the preburner station was 
also added to ailow mixing of the preburner 
exhaust gasses with the oncoming air so as not to 
suffocate the main burner jets of oxygen. A 
schematic showing the scramjet concept is 
presented in Fig. 2. 

Preliminary HAVOC design results at a Mach 
number of 8 and a dynamic pressure of 1500psf 
for the waverider HRV with two ramps and with 
both shock on shoulder and on cowl lip indicated 
that a contraction ratio of roughiy 14 was 
achievable. Given that the spontaneous ignition 
temperature of an ethylene/oxygen mixture is 
roughiy 780°K[17], a guideline of 1000°K was 
used as the design temperature at the 
mainbumer station to insure that combustion was 
indeed present. The 1-0 cycle code was run 
with the LOX augmentation prebuming option to 
compute parametrically, the required fuel and 
oxygen flow to achieve an equivalent mixed 1-D 
temperature at the main fuel injector station equal 
to the auto-ignition value. For an engine with an 
equivalence ratio of 1 (stoichiometric) this 
resulted in roughly 2.5% of the fuel being 
directed to the preburner that was then burned 
stoichiometricafy with onboard LOX. The flow of 
LOX was then reduced to the preburner 
resulting in a fuel-rich preburner. This reduced 
the required amount of onboard LOX with only a 
minimal effect on temperature distribution prior to 
main fuel injection. This is because the fuel-rich 
prebumer exhaust gasses continue to react with 
the air after injection. A heat balance on the 
vehicle and engine was used to compute the fuel 
total temperature. The resulting operating 
conditions of this engine are given in Table 1. 


Table 1 . Engine Operating Conditions 



Inlet 

Prebumer 

Main 

Gas 

Air 

Products 
and C 2 H 4 

C 2 H 4 

Mach No. 

3.83 

1.1 

2.2 

naaiH 

833 

3395 

801 

P(atm) 

0.86 

13.39 

2.49 

Area(in) 

2.0000 

0.005 

0.166 

Anale 

- 

90° 

30° 


The area in the above table is based on a one 
inch width section. 


An important consideration in the design of a 
scramjet is the penetration distance of the 


injectors. Reference 1 8 gives a relationship for 
the jet penetration distance as a function of the 
jet and free stream Mach numbers, the 
momentum ratio of the two streams, the angle of 
injection and the jet diameter. The prebumer 
injection was designed to only penetrate 
through the boundary layer, whereas the main 
injection was designed to reach well into the air 
stream for better mixing. This resuited in a 
guideline of h=0.5cm for the prebumer and 
h=2.0cm on the main burner. Because of the 
large amounts of fuel being injected through the 
main injectors at stoichiometric conditions, the 
main fuel injectors were designed with a 
streamwise length-to- width aspect ratio of 5 and 
an angle of injection of 30 degrees to help 
reduce blockage. The injectors were laterally 
spaced one inch apart on both the top and 
bottom of the scramjet. The top injectors were 
then offset one-half inch to a produce a 
staggered injection for the purpose of increasing 
jet penetration and to avoid the additional losses 
of impacting jets. The prebumers were aligned 
with the main fuel injectors to ensure that the hot 
prebumer gasses fell in near vicinity of the main 
fuel. 

A three-dimensional CFD analysis of this engine 
was performed. Because of the periodicity of the 
engine in the absence of any side walls, only the 
flow between the centerline of the top injectors 
and the centerline of the bottom injectors was 
actually solved. The grid for this computation was 
generated algebraically and contained 119 cells 
in the axial direction, 60 cells from top to bottom, 
and 16 cells laterally. The grid spacing on the 
walls was set to 4.0 x 10’ 5 m. All of the injector 
exits were modeled as rectangles containing 8 
ceils in each of the axial and lateral directions. 
This three-dimensional computation required 
1 850 iterations leading to 1 02 hours on a Cray- 
YMP. The computation was halted after no 
piotable difference was seen with further 
iteration. Throughout the CFD simulation 
process, the ingestion of a thick boundary layer 
formed on the forebody of the hypersonic 
vehicle by the scramjet was neglected. 

For the current engine design, all of the injectors 
(including prebumer and main) were designed to 
be supersonic in the boundary-normal direction. 
This simplifies the CFD boundary condition 
procedure since merely specification of the 
injector variables is required for a supersonic 
inflow boundary condition. This type of injector is 
also quite practical in a scramjet since injector 
pressures are typically high enough to choke the 
injector flow. The supersonic inflow boundary 
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condition was imposed oniy on those ceil faces 
that correspond to an injector exit. No-slip, 
viscous boundary conditions were imposed on 
the ceil faces adjacent to the injector exit. This 
led to the use of rectangular injectors to avoid 
further complication of the boundary condition 
procedure and the grid generation process. 

The results of the CFD analysis are shown in 
Figs. 3-13. The pressure contours on both the 
symmetry plane containing the top injector and 
the one containing the lower injector are shown 
in Fig. 3. The pressures are smoothly varying 
except in the vicinity of the injectors. Shocks 
emanating from both the preburner and main 
fuel injectors traverse the height and width of the 
computational space. These shocks can 
significantly affect the efficiency of the engine, 
and any further refinement of this design wiil 
address the losses associated with these 
structures. The temperature contours (Fig. 4) 
show the injector penetration and the 
temperature rise caused by combustion. This 
figure indicates that the penetration of the 
combustion region is more than half the height 
of the scramjet. A significant degree of 
penetration is present, without adverse effects 
such as jet-jet interactions and jet-wall 
interactions. 

The velocity vectors of Figs. 5 and 6 show 
details of the injector flow fields. The velocity 
vectors in the immediate vicinity of the prebumer 
exit exhibit a spreading behavior that is typical of 
an under-expanded jet. This phenomenon is 
also present in the main injector region but is 
less visible because of the inclination of the 
vectors. A separation region is present in the 
prebumer injector region that reaches seven jet 
diameters upstream of the injector . This 
phenomenon is absent near the main fuel 
injectors because of the reduced angle of 
injection. 

Figures 7, 8 and 9 contain crossflow contour 
plots of temperature, fuel and water 
respectively, at various axial locations. The axial 
locations correspond to the following: 1 ) just aft 
of the prebumer station, 2) just upstream of the 
main fuel injection, 3) the back of the main 
injector station, 4) within the combustion 
chamber, and 5) the combustor exit. The exact 
axial locations are indicated on the plots. 

The temperature contours of Fig. 7 clearly show 
the mechanism that is studied in this paper. Fig. 
7(a) shows the hot preburner gasses that 
emerge from the preburner injector ports. 


These gasses mix and react with the oncoming 
air stream but still contain a very hot core just 
before main fuel injection (Fig. 7(b)). This hot 
core, failing just above the main fuel injection, 
serves as a "pilot light" for main fuel injectors 
causing combustion of the main fuel to 
instantaneously occur (Fig. 7(c)). The main fuel 
injectors were designed to produce a significant 
amount of penetration without traversing the 
entire height of the scramjet. This was 
accomplished and is clearly shown in the 
combustion chamber temperature contours 
(Figs. 7(d&e)). These figures indicate that the 
concept of prebuming does indeed accomplish 
the task of maintaining combustion at the main 
fuel injection station and that an injector can be 
designed to provide significant flow path 
penetration without unstarting the engine. 

The effect of the upper surface corner on 
combustion is shown in Figs. 8 and 9. An 
expansion wave emanating from this corner 
causes the density and temperature to decrease 
having an adverse effect on the rate of 
combustion and mixing. This wave affects the 
upper gasses before reaching the lower gasses. 
Therefore the expansion has a greater effect on 
the upper gasses. For this reason, there are 
more unbumed and unmixed gasses present in 
the upper region of the scramjet. 

The remaining figures present a comparison of 
the CFD results with those of a 1-D cycle 
analysis. The combustor efficiency computed 
by the CFD solution was implemented in the 
cycle code since no other schedule was 
available for this engine design. This was 
accomplished by curve fitting the average fuel 
fraction schedule (Fig. 10) and using this 
schedule in the 1-D cycle code. Both the CFD 
predicted schedule and the curve fit are shown 
in Fig. 10. Also shown on Fig. 10 is the 
predicted amount of carbon monoxide from both 
the CFD and 1-D cycle analysis. The CFD 
analysis predicts a higher amount than the 1-D 
cycle code. This is caused by the difference in 
the equilibrium mechanisms of the two codes 
and suggests an improvement to the simplified 
kinetics employed in the CFD solver for scramjet 
computations. Comparisons of the average 
temperature, the momentum-averaged pressure 
and the mass-averaged velocity from the CFD 
analysis with the results of the cycle code show 
general agreement (Figs. 11-13). The 
discrepancies can be attributed to the improved 
ability of the CFD method to account for detail 
and the difference in the equilibrium 
mechanisms. Finally, a sensitivity analysis using 
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the cycle code indicates that, at Mach 8, a 1% 
change in overall combustion efficiency 
represents approximately 1% change in cowi-to- 
taii Isp and 0.8% change in axial thrust 
coefficient. 

Hypersonic Waverid^r RQS9arQ tLV9hislg-DfisL<aD 
and Analysis 

Unlike acceierator-type missions (e.g., SSTO) 
where the mass capture characteristics of the 
vehicle are most important, cruise configurations 
place an emphasis on the aerodynamic 
performance of the design. Hence vehicle lift-to- 
drag ratio (L/D) or the product of L/D and engine 
Isp (a parameter proportional to Brequet range 
factor) becomes more important. A waverider 
configuration, with high hypersonic L/D potential, 
was selected as the baseline configuration for the 
present study of the HRV. The 
configuration/engine installation was numerically 
optimized to maximize (L/D) x Isp, using 
forebody shape, ramp angles and cowl position 
as optimization parameters. For this analysis, the 
forebody is a Mach 8 waverider configuration and 
the ramp angles and cowl position were designed 
to produce first and second ramp shock-on-lip 
and cowl shock-on-shoulder. The engine 
geometry and operating parameters were held 
fixed in the optimization procedure. The 
numerical optimization was performed using the 
HAVOC cycle code. 

The design parameters and constraint functions 
used for optimization are listed below. There 
were 15 design parameters used in the 
optimization process. Six parameters defined the 
waverider generating surface (hence the vehicle 
shape); eight parameters defined the ramps, 
cowl, nozzle geometries, and leading-edge 
radius; the last design parameter was free stream 
dynamic pressure. Eight constraints were used 
in the optimization process, including engine 
throttle greater than or equal to the required 
cooling equivalence ratio, leading edge 
equilibrium radiation temperature limited to 
3200° F, wing tip closure angle limited to 10 
degrees, and the vehicle structural thickness limit 
at vehicle aft end, accounting for nozzle 
integration. The volume was constrained to 1 .7% 
of vehicle length cubed. Finally, the vehicle 
width to length ratio was less than or equal to 
0.75. 

For the Mach 8 design, a generating shock angle 
of 12 degrees was arbitrarily selected. The 
waverider shape optimization process involved 
planform shape changes to sweep the leading 


edge in order to alleviate high heating rates at the 
higher free stream dynamic pressures, traded off 
against leading edge radius and associated 
leading edge bluntness drag. The design 
optimization process of the inlet (ramp positions 
and ramp angles) resulted in the two ramp shocks 
converging on the cowl lip (shock-on-lip) then 
reflecting on the shoulder of the combustor 
entrance (shock-on-shoulder). The geometric 
contraction ratio was approximately 14, with a 
resulting pressure at the combustor iniet of 
about one atmosphere. Preliminary performance 
estimates indicated a required vehicle length of 
23 ft, with an overall vehicle body density of 20 
lbs/ft**3. A lift-to-drag ratio of 4.3 was achieved, 
with a cowl-to-tail Isp of 746 seconds at an 
assumed combustor efficiency of 95%. and at an 
equivalence ratio of 1.0, with a free stream 
dynamic pressure of 900 psf, resulting in an initial 
cruise altitude of 92,500 ft. An engine width of 
roughly 0.762m produced a net thrust that is 
equal to the net drag of the vehicle at the design 
point. For the results presented below, the 
engine combustion efficiency of 55% was taken 
from the scramjet engine CFD results. 

Nose-to-tail analyses utilizing both CFD and cycle 
codes were performed. The geometries for both 
analyses, however, were slightly different. The 
ramp sidewalls in the cycle analysis were 
assumed to have zero thickness and were 
aligned with the forebody flow. They therefore 
resulted in zero pressure drag with only a nominal 
addition to the viscous drag. The CFD geometry, 
on the other hand, contained ramp sidewalls with 
an included angle of 21 degrees Because the 
sidewall attaches to the cowl lip and because the 
cowl protruded far below the lower surface of the 
waverider, the sidewall contained a significant 
amount of additional volume. This added volume 
was exposed to the forebody flow and produced 
additional drag to the vehicle. 

The grid generation for the CFD nose-to-tail 
analysis is described below. The grid generation 
tools used for the external grids are the S3D and 
HYPGEN-UI codes mentioned previously. The 
external portion was divided into three axial 
sections; nose-to-inlet, inlet-to-exit, and exit-to- 
end sections. Grid generation was performed 
separately on each of these portions. For the 
surface grid, points were clustered on the 
waverider leading edge, on the sidewall leading 
edge and on all of the other convex comers of 
the vehicle. The lower surface grid is shown in 
Rg. 1 4. The outer boundary of the volume grid 
was placed well beyond the anticipated bow 
shock to a distance of roughly three meters from 
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the vehicle surface. The spacing of the first point 
from the surface was set to 1.0x1 0' 5 m. The grid 
dimensions for the three externa! grids measured 
102x70x144, 60x70x129 and 12x70x130 from 
nose to tail. The values of these dimensions 
correspond to the number of grid points in the 
streamwise direction, circumferential direction, 
and radial direction. An internal nozzle grid 
measuring 3x79x79 was generated algebraically 
for the CFD solution. 

The space-marching scheme, STUFF, was 
employed to obtain the nozzle and external CFD 
results for the waverider HRV. The space- 
marching solution began at the nose of the 
vehicle by setting the dependent variables in all 
of the cells to be free stream. The solution 
progressed by marching downstream through 
each of the grids until the aft end of the vehicle 
was reached. A fully conservative patching 
scheme [19] was employed to transfer the 
solution from one grid to the next. The frozen 
chemistry option was employed for the CFD 
results including the nozzle portion. This was 
sufficient since combustion is slowed 
considerably at the combustor exit. 

For the CFD nose-to-tai! results presented here, 
the HAVOC cycle code was used to provide the 
solution within the inlets and combustor. This 
was done to avoid the computational expense of 
obtaining time-marching solutions of both the 
outboard and inboard engines with 48 fuel 
injection ports in each. The CFD forebody 
results were averaged at the inlet face to provide 
one-dimensional inlet conditions for the cycle 
code. The presence of an oblique cowl shock 
was accounted for and provided the mechanism 
to turn the flow parallel to the engine cowl. The 
CFD-predicted heat release schedule presented 
above was used in the one-dimensional cycle 
analysis of the combustor. 

The two-dimensional, shock/weak wave results of 
the HAVOC analysis are plotted in Fig. 15. This 
figure shows the shock, contact, and expansion 
wave arrangement for the current integrated 
design on the vehicle symmetry plane. The 
intent of the forebody/ ramp design is clearly 
shown in this figure. Both of the ramp shocks 
impinge on the cowl leading edge and the cowl 
shock impinges on the shoulder. The bow 
shock, however, lies outside of the inlet. This is 
because at Mach 8, with a fixed throat height, 
having the bow shock on the cowl lip can 
produce shock-on-ramp with possible boundary 
layer separation and engine unstart. The nozzle 
geometry was restricted to two planar sections 


and the corresponding expansion and shock are 
shown in Rg. 15. 

Rgure 16 is a line plot of the surface pressure on 
the keel line of the vehicle. Both the CFD and 
cycle results are plotted. Agreement is excellent 
except in the vicinity of inviscid/viscous 
interactions. Comparison of the pressure on the 
waverider forebody is within 2 percent except at 
the leading edge where CFD predicts a pressure 
spike resulting from boundary layer 
displacement. The ramp pressures also agree 
quite well except for the asymptotic behavior of 
the CFD results that is typical of shock/boundary 
layer interactions. Another reason for this 
discrepancy can be attributed to the three- 
dimensional intersection of the first ramp with the 
curved waverider forebody. The intersection line 
bends downstream from the vehicle symmetry 
plane. This three-dimensional effect reduces the 
pressure near the intersection region. Rgure 16 
also shows very good agreement on the outer 
surface of the cowl and on the nozzle surface. 
Comparison of surface pressure on the upper 
surface showed that the CFD predicted 
pressures were 3-10% higher than the 
engineering predictions. . This is explained by 
viscous/inviscid interactions. 

The pressure contours in Rg. 17 show the keel 
line shock and expansion waves that are 
predicted with CFD. The location and strength of 
these waves are in very good agreement with 
those predicted by the cycle code (Rg. 15). The 
CFD results, however, predict that the ramp 
shocks lie slightly outside of the inlet instead of 
on the cowl lip. This difference is attributed to 
shock/boundary layer interactions and to three- 
dimensional effects. 

Rgure 18 shows the CFD predicted density 
contours on crossflow planes at various axial 
stations. These stations correspond to the 
following locations: a) on the second ramp, b) at 
the inlet, c) at the middle of the cowl, d) at the 
engine exit and e) at the vehicle end. The two 
ramp shocks are visible in the pressure contours 
of Fig. 18a and clearly exhibit a 3-dimensional 
behavior. The bow shock is attached to the 
leading edge of the waverider and produces a 
very clean flow field on the forebody absent of 
the ramp geometry. 

Figure 18 shows that the flow into the inlet is 
fairly clean except for a weak shock that emanates 
from the ramp sidewall and interacts with the 
forebody boundary layer. Also evident is a 
thickening of the ramp boundary layer thickness 
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nearer the ramp sidewall. This is caused by the 
tendency of the flow to spill off the first ramp 
resulting in a sideways velocity component. This 
effect washes the boundary layer away from the 
symmetry plane of the vehicle. It is then stopped 
before the inlet station by the sidewalls on the 
second ramp. Another feature that should be 
noted in Fig. 1 8b is the shock emanating from the 
outer surface of the ramp sidewall. This feature is 
absent in the cycle analysis since the sidewalls 
were assumed to have zero thickness. The 
resulting high pressure region adds considerably 
to the net drag of the vehicle since the outer 
sidewall has an area component in the 
streamwise direction. This effect resulted in a 
negative net-thrust predicted by the nose-to-tail 
analysis, and therefore, creates a need to reduce 
the thickness of the sidewalls. At the cowl mid- 
point (Fig. 18c), the ramp shocks have clearly 
merged with the bow shock an expansion fan 
emanating from the lower surface of the cowl is 
evident. 

Figure 18(d-e) depict the CFD predicted flow field 
on the aft portion of the vehicle. This flow field 
region is very complex due to the presence of 
the engine, although a significant portion of the 
wave rider flow field remains unaffected. The bow 
shock remains attached to the leading edge for 
the entire length of the vehicle. This feature is 
desirable since any spillage of the high pressure 
gasses onto the upper surface would reduce the 
performance of the vehicle. Two flow field 
features that are clearly visible in the nozzle 
portion of the flow are the initial nozzle expansion 
and the shock caused by the second nozzle 
plane. These features appear to be nearly two- 
dimensional near the centerline, leading to the 
good agreement with the cycle analysis (Fig. 16). 

Figure 19 shows the pressure contours on the 
lower surface of the vehicle. The impact on the 
surface pressure by the presence of the 
integrated scramjet is visible in this figure. The 
lower surface is composed of a large region of 
undisturbed flow with a fairly constant pressure. 
The remaining surface is exposed to the 
numerous shocks and expansion fans that 
originate from the integrated propulsion system. 
These include the ramp shocks, the sidewall 
shocks, the exterior cowl expansions, the engine 
shroud expansion, and the numerous nozzle 
waves. 

Of particular interest in the current design are the 
forebody ramp system and the resulting flow 
field. Due to the curvature of the waverider 
forebody, an oblique intersection results 


between it and the first ramp. This intersection is 
visible in the pressure contours of Fig. 1 9. The 
strong pressure gradients that are present in this 
vicinity cause the flow to be diverted away from 
the vehicle symmetry plane beginning at the 
origin of the first ramp. This effect tends to 
reduce the amount of mass provided to the inlet. 
An intersection that is less oblique would reduce 
the spillage since the pressure gradient would be 
aligned more with the flow direction. This can be 
accomplished by reducing the curvature of the 
waverider forebody at the symmetry plane. In the 
current design, this spillage was effectively halted 
at the second ramp by using sidewalls. This is 
shown in the surface streamlines of Fig. 20. 

CONCLUSIONS 

Two methods were demonstrated for the 
conceptual design and analysis of a hypersonic 
waverider research vehicle: an engineering 
analysis code with a simplified nose-to-tai! flow 
field analysis capability and a complete 3-D CFD 
code with a hydrocarbon/air capability. The 
methods have been shown to produce good 
agreement for waverider forebodies except in 
regions of strong viscous/inviscid interaction. 
The coupling of the cycle code with CFD to 
compute a single flow field was also 
demonstrated. Even though the cycle code has 
significantly less ability to predict detail, this 
coupling proved to be useful because it reduced 
computer requirements. 

The nose-to-tail analysis of the waverider HRV 
clearly shown the benefits of the current design 
and revealed areas for improvement. The 
waverider forebody combines with the current 
ramp system to provide a uniform flow to the inlet 
of the scramjet. Further, the forebody shock 
remains attached to the leading edge for the 
entire length of the vehicle. This avoidsspillage 
of the high pressure air onto the upper surface 
that can significantly reduce lift. The analysis has 
also shown that the ramp sidewalls need to be 
thinner since they protrude far into the forebody 
flow field and cause unnecessary drag to the 
vehicle. 

The analysis of the liquid oxygen-augmented 
prebuming hydrocarbon scramjet indicates that 
the concept does indeed produce combustion of 
the main fuel within the scramjet engine. The 
preburning process provides a sufficiently 
elevated temperature flow into the main fuel 
injector region to support immediate combustion 
of the gaseous ethylene fuel. However, 
because of the significant amount of unbumed 
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fuel at the combustor exit, there is a need for 
better mixing efficiency within the combustor. 
For the current fuel injector configuration, 
improved engine cycie performance could also 
be achieved at a lower overall engine 
equivalence ratio, limited by engine cooling 
requirements. Finally, there is currently a void of 
high-quality, CFD validation-type data for high- 
speed hydrocarbon combustion, leaving an 
uncertainty in any predicted results. 
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Fig. 1 Hypersonic vehicle concept with an 
integrated hydrocarbon scramjet 
engine. 



Fig. 2 Schematic diagram of the 3-D scramjet 
engine geometry 



Fig. 4 Temperature contours at a) symmetry 
plane of top injectors, b) symmetry 
plane of bottom injectors. 



Fig. 5 Velocity vectors in vicinity of prebumer 
injector 



Ftg. 3 Pressure contours at a) symmetry plane 
of top injectors, b) symmetry plane of 
bottom injectors. 



Fig. 6 Velocity vectors in vicinity of main 
burner injector 
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Figure 7 Temperature Contours at indicated Axial Locations 
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Figure 8 C 2 H 4 Mass Fractions Contour at Indicated Axial Locations 
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H 2 O Mass Fraction Contours at Indicated Axial Locations 
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Fig. io C2H4 and CO mass concentrations as a 
function of axial location 
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Fig. 13 Axial velocity as a function of axial 
location 
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Fig. 11 Temperature as a function of axial 
location 


Fig. 14 Lower surface grid for the current 
waverider HRV design. 



— TUFF 
o HAVOC 


Fig. 12 Pressure as a function of axial location 


Fig. 15 HAVOC predicted flow field on the 
lower symmetry plane. 
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Fig. 16 Keel line pressure comparison 
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Fig. 17 CFD predicted pressure contours on 
vehicle symmetry plane. 
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Fig. 19 CFD predicted surface pressure plot. 
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Fig. 20 CFD predicted surface streamlines on 
the forebody ramp system. 


Fig. 18 CFD predicted density contours at 
various axial locations. 
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ABSTRACT 

Computations are done on flow configurations that 
resemble the reaction zone in the scramjet combus- 
tor flows. Compressible, reacting, turbulent flow so- 
lutions are obtained. A two equation (I-c) model 
with compressibility correction is used to calculate the 
flow field. A finite rate (8-species, 13-reaction steps) 
chemistry model for hydrogen-air combustion has been 
used. Computations are carried out using the Navier- 
Stokes solver TUFF. Predictions are compared with 
available experimental data and also those obtained 
by using the code UPS. 
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stress tensor 
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Subscripts 


A,b 

coefficients in Arrhenius rate equation 


turbulence model constants 
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total internal energy 

fn 

mass fraction of species n 

H 

total enthalpy 

h 

static enthalpy 
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turbulent kinetic energy 


forward and backward reaction rate constants 
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number of reaction steps 

M n 

Molecular weight of species n 

M t 

Turbulent Mach number 
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number of chemical species 

Pr,Pr t 

laminar and turbulent Prandtl numbers 
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pressure 

Sc,Sc t 

laminar and turbulent Schmidt numbers 
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temperature 

T a 

Activation Temperature 
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time 
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velocity vector 

td n 

production rate of species n 
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streamwise coordinate 
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INTRODUCTION 

Hypersonic travel requires propulsion systems which 
are different from the conventional ones used in most 
of the modem aircraft. The supersonic combustion 
ramjet (scramjet) is a system considered to be suit- 
able for high speed applications. There has been a 
tremendous amount of activity in the area of scramjet 
research in recent years ( [1] - [11]). Some of the re- 
lated topics include inlet configuration, mixing layers, 
mixing enhancement, combustor configuration, finite 
rate chemistry models and chemical kinetics. The fuel 
used in the scramjet varies depending upon the ap- 
plication. For example, hypersonic waveriders using 
hydrocarbon fuels have been designed [12] for applica- 
tions in the moderate hypersonic speed regimes. For 
Mach numbers of the order of 15 and above, hydro- 
gen is generally considered to be the fuel of choice. 
In the present work, hydrogen is the fuel used in the 
computations. 

The present work represents a computational ef- 


1 



fort in establishing a solution procedure for hypersonic 
propulsion applications. The entire task of establish- 
ing the solutions procedure must then be divided into 
smaller tasks dealing with subsets such as turbulence 
modelling, chemical kinetics, geometry etc. One such 
task is the topic for the present study. Here, the rele- 
vant flow features of the combustor, namely the mixing 
and chemical reaction between the fuel and oxidizer 
streams, is addressed. The flow field in a scramjet is 
complex. It is turbulent and compressible involving 
high heat release. The solution procedure should ad- 
dress all aspects of the flow field adequately. It should 
be capable of accurately modelling the turbulent field, 
taking into account the effects of compressibility, and 
addressing the changes associated with heat release. 
Also, the interactions between the distinct physical 
aspects of the flow such as the effect of heat release 
on turbulence, the interaction between turbulence and 
chemistry etc. must be properly addressed. Signifi- 
cant progress has been made in addressing these areas 
via accurate and realistic modelling in recent years [6]. 

Remarkable advances have been made in the area 
of turbulence modelling, accounting for a variety of 
factors that affect the flow field. Compressibility cor- 
rection models to account for the effects of compress- 
ibility, near-wall turbulence models to deal with the 
transition from fully turbulent to zero turbulence, vis- 
cous dominated flow field near no-slip boundaries and 
modifications to models to account for flow curvature 
are some examples. A wide variety of turbulence mod- 
els, including algebraic (zero-equation), one-equation, 
two-equation, Reynolds stress and large eddy simula- 
tion models, are available (for example references [13] - 
[15]) depending upon the sophistication and accuracy 
desired and the limits imposed by numerical solution 
procedures. 

Thermodynamic and chemical kinetic models [16] 
applicable to the scramjet flows have been undergoing 
continuous improvements in recent years. Accurate 
modelling of thermodynamic variables as functions of 
temperature which are valid over a wide range of tem- 
peratures is an example. In flows such as the one 
associated with the scramjet, the time scales associ- 
ated with fluid dynamics and chemical reaction (not to 
mention the turbulence scales) require that the com- 
bustion process be modelled via a finite rate chem- 
istry mechanism. Such a mechanism should account 
not only for the major species (reactants and prod- 
ucts) involved in the chemical reactions but also the 
intermediate transient ones which play a vital role in 
the reaction progress process. Accurate models for the 
chemical reactions in the scramjet combustor, thus, is 
a crucial aspect of the solution procedure. 

The design of the combustor is strongly dependent 


upon factors such as the mixing between fuel and ox- 
idizer streams, presence of shocks in the flow field, 
boundary layer effects, flow separation, extent of chem- 
ical reaction within the combustor and so on. The nu- 
merical solution procedure should have the capability 
of addressing all of these factors while maintaining the 
required accuracy and robustness. There is a glut of 
useful numerical solvers applicabtf |j|jt a wide variety 
of flows including all speed regimes.. Computational 
algorithms which are fast and accurate are being im- 
proved everyday. 

Even though there are a wide variety of sophisti- 
cated and physically accurate thermodynamic, chem- 
ical kinetic and turbulence models available it is not 
always possible to use the most accurate and elabo- 
rate versions in a numerical simulation due to the lim- 
itations imposed by computer memory requirements, 
computational economy, ease of use and adaptability 
to practical problems. Solutions often are required, 
especially in ihe engineering industry which is the end 
user for such solvers, in a short time using comput- 
ers that may not be the fastest available. As a result, 
compromises must be struck between physical accu- 
racy and computational feasibility and it is this aspect 
which differentiates between various solvers that exist 
today. 

In the present study, an attempt is made to es- 
tablish a solution procedure for scramjet combustor 
flow predictions from the perspective of the discus- 
sion above. The models chosen to represent the tur- 
bulent and chemistry fields reflect the compromise be- 
tween physical accuracy and computational economy 
mentioned above. The code chosen for the computa- 
tions is the TUFF [17] code and the solutions are com- 
pared with those obtained with the UPS [18, 20] code. 
The turbulence model chosen is the two-equation k — c 
turbulence model with low Reynolds number modifi- 
cations [13]. However, the Baldwin- Lomax algebraic 
model is also available as an option. The compress- 
ibility effects are included via the compressibility cor- 
rection model proposed by Zeman [19]. The fuel used 
is hydrogen although the numerical solver can easily 
be modified for hydrocarbon fuels. A 9-species, 20- 
reaction steps chemical kinetics model for hydrogen- 
air combustion [16] is available. For the computations 
presented in this report, an abbreviated version (8- 
species, 13-steps) of this model has been used. 

Mixing plays a major role in high speed combus- 
tor flows. The reaction zone is mainly confined to 
mixing layers that exist between fuel and oxidizer 
streams. Two flow configurations are chosen for the 
study. The first is the well known Burrows-Kurkov 
experiment [21] in which hydrogen and vitiated air 
streams (two-dimensional) mix and react. The second 
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case is that of an axisymmetric configuration (22, 23] 
where two coaxiai jets (fuel and oxidizer) mix and re- 
act. Experimental data from compressible, reacting 
mixing layers is still scarce which hinders the valida- 
tion of the calculation procedure. The available data 
from the above two experiments are used to compare 
with the predictions. The governing and secondary 
equations used in the computations have all been de- 
scribed in detail in the references cited above. Only an 
abreviated equation set will be given in the present pam- 
per. The computations were performed on the super- 
computers of NAS and NASA Ames Research Center 
(090). 


GOVERNING EQUATIONS 

The equations used for computations are described 
in detail in references [6, 7, 10] and [24]. Only the 

forms of the modelled equations used in the present 

study are given here. Density-weighted averaging is 
used to derive the mean flow equations from the in- 
stantaneous coservation equations. The dependent 
variables, with the exception of density and pressure, 
are written as 

<f> = <f> + <j>" (1) 

where the is the fluctuating component of the vari- 
able under consideration and its Favre-mean 4> is de- 
fined as 

1> = P 4 (2) 

P 

In this equation, the overbar indicates conventional 
time- averaging. Density and pressure axe split in the 
conventional sense as, 


and 


e 


<9u" du" 

P v st; Tt; 


( 8 ) 


Boussinesq approximation is used to obtain closure 
of the averaged equations. Hetq the Reynolds stress 
tensor is written as, 


Sij 


. - pu'/u'j' = fi t Si, - ^ pkSij 

dUj dU) 2 dUl 

dxj dxi 3 dxt ,J 


(9) 


where fit is the turbulent/eddy viscosity defined as 


fj t = Cpp 


k 7 


( 10 ) 


with C^=0.09. 

The modelled momentum equation, then, is, 

dpVj dpJTJTj dp 2 dpk 

dt + dxj dxi 3 dxj lJ 

- os + *)£■)] (ii) 

The effects of compressibility are included via the 
model proposed by Zeman [19]. Here, the compress- 
ible dissipation terms axe expressed as functions of the 
turbulent kinetic energy dissipation rate and the local 
turbulent Mach number. The compressibility effects 
axe represented by a component of the dissipation rate 
(e c ) given as 


p — p + p and p = p + p' (3) 

The averaged continuity and momentum equations axe 

diUi 


2 + 

dt 


dxi 


= 0 
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dpUi , 
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dpu'lxi'! 

dt 

dxj 

dxi 
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with repeated indices indicating summation. 

In the two-equation turbulence model, the two tur- 
bulence variables axe the turbulent kinetic energy (£) 
and the dissipation rate (c) [13] defined as 


k = 


_ pu'fu'l 


2 P 


CO 


e« = Kct 
K e = tjF(Mt) 
2k 


Mt ~ 


|2 i 


( 12 ) 


where a is the local speed of sound and F(M t ) is a 
function of the local turbulent Mach number {Mt). 
F{Mt) is given by 


F(M,) = l-«p[-( M ' 06 M, ! ) 2 l, M,>M„ 

— 0, M t < Mt # 

with M t ^=0.1 and q=0.75. The modelled turbulent 
kinetic energy equation is (6, 13] 


dpk dpkUj 
dt + dxj 


= P k -pt(l + K c ) 

d Pt . dk 

+ x — l(p H )j— 
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(13) 
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m, 


,v 

E 
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i* 



1 = 1 


The modelled i-equation used (no compressibility 
corrections) in the present analysis [13] is given below. 


dpi dptUj 
di + dxj 


= (Ci Pi ~ C?pe) - 

d ■ _ fj t .dc. 


(15) 


where, i/'j and u" t axe the number of molecules of the 
scalar 5 involved in the /-th reaction step in the forward 
and backward directions, respectively. The forward 
and backward rate-constants of the reaction l are given 
by k/i and jbjj respectively. 

h„ = AT*' «p( -^] (21) 


where Pi is the production term in the turbulent 
kinetic energy equation. The model constants used 
in the analysis are Ci = 1.44, C 2 = l-92, <ri = 1.0, 

<r f = 1.3, Pr=0.72, Pr, = 1.0, 5c=0.22 and 5c, = 1.0. 

The mass-averaged total energy can be written in 
terms of the total enthalpy as 


E ~ H (16) 

P 

The correlations between the fluctuating velocity and 
the scalar fluctuations are modelled using a gradient- 
diffusion hypothesis. A typical model is of the form 


- puP<p" - 


Pt 
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(17) 


where <j± is a coefficient which, normally, is a constant. 
For 4> — f n (n represents the species) , cr # = Sc t , and 
for the static enthalpy, {<f> = h), <x $ = Pr t . Using the 
above definition, and omitting the body force contri- 
bution, the time- averaged and modelled energy equa- 
tion [6] is 


dpE 


+ = /-(W 

OXj OZj 


+ dTj Er + Pr,'dx, 




where < 7 * comes from the turbulent kinetic energy 
equation. The modeled species continuity equation is 


dpf n dpfnUj 

di dxj 


!; + &>§£i< 19 > 


The modelled form of the mean species production 
rate due to chemical reaction (tu„) is given, for a finite- 
rate system involving L reaction steps and N species, 
in the following general form: 


w n = T>”, ~ U nl) x 


(20) 


?=i 


j=i *— t 


where A;, b f and T ot are numerical constants specific 
to the given reaction step /. k n is determined from the 
equilibrium constant for the /-th reaction step and kj j. 

Solution of the modeled equations 

The equations are discretized and integrated in 
space and time to obtain steady state solutions using 
the finite- volume based numerical solver TUFF [17]. 
The TUFF code contains many desirable features for 
the computation of three-dimensional, hypersonic flow 
fields. It has non-equilibrium, equilibrium and perfect 
gas capabilities along with an incompressible option. 
It employs a finite-volume philosophy to ensure that 
the schemes are fully conservative. The upwind in- 
viscid fluxes are obtained by employing a new tem- 
poral Riemann solver that fully accounts for the gas 
model used. This property allows the flow field dis- 
continuities such as shocks and contact surfaces to be 
captured by the numerical scheme without smearing. 
Total Variation Diminishing (TVD) techniques are in- 
cluded to allow extension of the schemes to higher or- 
ders of accuracy without introducing spurious oscilla- 
tions. The schemes employ a strong coupling between 
the fluid dynamic and species conservation equations 
and are made fully implicit to eliminate the step-size 
restriction of explicit schemes. This \i necessary since 
step-sizes in a viscous, chemically reacting calculation 
can be excessively small for an explicit scheme, and 
the resulting computer times prohibitively large. A 
fully conservative zonal scheme has been implemented 
to allow solutions of very complex problems. The 
schemes are made implicit by fully linearizing all of 
the fluxes and source terms and by employing a mod- 
ified Newton iteration to eliminate any linearization 
and approximate factorization errors that might oc- 
cur. Approximate factorization is then employed to 
avoid solving many enormous banded matrices. As 
mentioned before, the options for turbulence models 
include both zero and two equation models (both k — e 
and k — w). For more details abqut the solution pro- 
cedure the reader is directed to the reference cited 
above [17]. 
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RESULTS AND DISCUSSION 

Two reacting flow configurations have been chosen 
for the present study. As mentioned before, the Navier 
Stokes solver TUFF has been used for the computa- 
tions. The first one is the case of coaxial jets (22, 23] 
where a hydrogen jet flows (inner jet) coaxially with 
an outer vitiated air (mass fractions: oxygen=0.246, 
water=0.209 and nitrogen=0.545) jet. A schematic of 
the flow problem is given in Figure 1. The two streams 
are, air (£7=1380 m/sec, T=1180 K with p=107000 
N/m2) and hydrogen (£7=1774 m/sec, T =545 K with 
P = H2000 N/m2)., The air stream is supersonic with a 
Mach number of 1.97 and the hydrogen stream Mach 
number is 1.00. The inlet mean velocity is assumed 
to have a step profile with the two jets having uni- 
form speeds at the specified values (no experimental 
data available). The velocity in the lip region of the 
inner jet tube wall (finite wall thickness) is assumed 
to be zero. The inlet temperature profile is derived 
based on the experimental data given for a location 
just downstream of it (shown later). The inlet species 
mass fraction distributions are also chosen based on 
the experimental data provided at the same down- 
stream location. A constant turbulence intensity level 
is used in the free stream for arriving at the initial 
distribution of turbulent kinetic energy and the dis- 
sipation rate. A 13-step, 8-species Hi — Air reaction 
model (Table 1) has been used for the finite-rate chem- 
istry system considered here. A 81 X 91 grid (81 points 
in flow direction, 91 points in the radial direction) was 
used for the calculations. The inner jet/tube diameter 
(D~Q. 00236 m) is used as a reference length. The to- 
tal length of the flow domain is equal to 43.1 D. The 
outer boundary (radial) of the flow domain is taken to 
be at y=17 D. A more detailed description of the flow 
parameters is given in Table 2. The region outside the 
limits of the air jet is assumed to be still air at a tem- 
perature of 273 K. The two-equation (it — <) turbulence 
model is used along with the finite rate Hi- Air chem- 
istry model mentioned above. In all the figures shown 
in this report, y refers to the radial distance measured 
from the axis of the coaxial system of jets. 

Figures 2-3 show the results of the computations. 
Figure 2 shows the computed and experimental distri- 
butions of species mole fractions. The figure is de- 
signed in a two-column format. The left side col- 
umn represents the inlet (first x-location) data and 
the right side column is the data at the exit plane 
(x/D=43.1 D ). As seen in these figures, the inlet 
data agreement between the computations and exper- 
iment is not perfect, especially around the jet edges, 
and this might affect the computed distributions at 
downstream locations. The comparison between pre- 
dictions and experiment at the downstream location 


{z/D 43.1 D) is good given the above mismatch be- 
tween the two data at the inlet. The development 
of the reaction zone after ignition is not predicted 
well. The experimental data indicates that the reac- 
tion zone (depicted by the water mole fraction distri- 
bution) spreads more quickly than the predictions in- 
dicate. The predictions show the reaction zone to be 
off-center whereas the experimental data shows the re- 
cation zone to be closer to the axis of symmetry. How- 
ever, there is very good qualitative agreement berween 
the data with the peak values of the reaction prod- 
ucts predicted very well. The flow domain was seen 
to have a wave-like structure as shown by the pre- 
dicted profiles. The worst agreement seems to be for 
the case of oxygen. However, when the initial pro- 
files of oxygen are compared one finds that there too 
is the worst agreement between computations and ex- 
periment which may be reason for the problem down- 
stream. Figure 3 shows the comparison of static tem- 
perature data. The agreement between predictions 
and experiment is good qualitatively displaying similar 
trends. The uncertainty associated with the accuracy 
of the experimental data is unknown. There are con- 
siderable differences between the data presented by the 
two references [22, 23], especially in the temperature 
profiles. Overall, there is good qualitative agreement 
between the predictions and experiment. 

The second test case considered is the Burrows- 
Kurkov experiment [21]. The flow configuration is two- 
dimensional. A schematic diagram of the configuration 
is given in figure 4. No-slip walls bound both the upper 
and lower regions (y=0 and y=y maar ). The lower wall 
is inclined to form an expansion surface. Hydrogen is 
injected along this surface into a vitiated air stream. 
The two streams mix and react downstream of the 
injection location (inlet). The hydrogen stream is in- 
jected at a velocity of 1216 m/3ec and a temperature 
of 254 K. The airstream comes in at a speed of 1764 
m/sec and a temperature of 1270 1^. Full details about 
the flow parameters and geometry are given in Table 3. 
In this case, the reference length used in the hydrogen 
jet width at inlet, h (A=0.0Q4 m). The models for tur- 
bulence and chemistry are identical to the ones used 
for the coaxial jet case. The grid size is 81 X 121 
(81 grid points in the axial (x) direction and 121 grid 
points in the transverse direction). The total length 
of the solution domain is 0.356 m (x/A=89). Avail- 
able inlet data have been used for the first-plane pro- 
files which improved the predictions remarkably over 
the solutions obtained with unif orm profiles. The so- 
lutions are compared with the available experimental 
data at this location (exit) in figures 5-7. The solu- 
tions carried out with the space marching PNS code 
UPS [20] using the Baldwin-Lomax turbulence model 
are also given for comparison. 
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Figure 5 shows the compariosn between the pre- 
dicted distributions of the species mole fractions and 
the corresponding experimental data. As seen in 
these figures, there is excellent agreement between the 
TUFF predictions and experiment. The predictions by 
the UPS code do not agree very well but still there is 
very good qualitative agreement with the experimental 
data. Figure 6 compares the predicted profiles of exit 
plane total temperature and Mach number with the 
experimental data. There is good qualitative agree- 
ment in the case of temperature and very good overall 
agreement in the case of the Mach number distribu- 
tion. Figure 7 shows the comparison between the pre- 
dictions and experiment of the lower wall (hydrogen 
jet side) pressure. Ignition causes the pressure rise in 
the profile. Ignition seems to be delayed in the case 
of the predictions accompanied by a more pronounced 
pressure rise. 

High speed reacting flows such as the two cases stud- 
ied here are complex inspite of their simple geome- 
tries. The interactions between the different aspects 
of the flow such as turbulence, chemical kinetics, heat 
release etc. are very difficult to understand and, to 
a large extent, impossible to model accurately. Mod- 
ern day experimental facilities still cannot make com- 
plete measurements in such flows. Only mean values 
of temperature, velocity, pressure, species concentra- 
tions etc. are available, if any. Even then, the un- 
certainties associated with the data force one to ac- 
cept them only with certain reservations. Given that, 
there is almost never a chance for perfect agreement 
between predictions and experiment in all the areas. 
While the advances made in measurement techniques 
improve every day, the fruits of these advancements 
(ie. accurate measurements) are not realized imme- 
diately. As a result, today's computations will have 
only old data for validation (in the present case, the 
best data is already four years old) purposes which is 
certainly the case here. Unless more accurate experi- 
mental data with minimum uncertainties are available, 
the best a computational effort can hope for, in terms 
of validation, is probably what is seen here. 

CONCLUSIONS 

Computation of the flow fields of two high-speed, 
turbulent, reacting flow configurations involving finite- 
rate chemical kinetics for hydrogen-air combustion 
have been carried out. A two-equation (k - c) turbu- 
lence model with compressibility corrections has been 
used. The predictions are compared with available 
experimental data. Good qualitative agreement is 
present between computations and experiment. More 
detailed experimental data is necessary. 
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Table 1. Hz — Air Reaction System 


No. 

Reaction 

1 

H + O 2 — O + OH 

2 

OH + H 2 ^ H 2 O 4- H 

3 

0 + Hz ^ OH + H 

4 

OH + OH — H 2 O + O 

5 

H + OH 4- M H 2 O+ M 

6 

H + H + M — H 2 +M 

7 

O + O + M ^ O 2 - h M 

8 

H + O+M — OH + M 

9 

H 4- O 2 4- M H O 2 4* 1W 

10 

OH 4- HO 2 — H 2 0 4- O 2 

11 

H 4- HO 2 — H 2 + O 2 

12 

H + HO 2 ^OH + OH 

13 

0 4* HO 2 ^ OH 4* O 2 

14 

HO 2 4- HO 2 — H 2 O 2 4- O 2 

15 

H 4- H 2 O 2 ^Hz + HO 2 

16 

OH 4- H 2 O 2 — H 7 0 4- HOz 

17 

H 4- H 2 O 2 — H 2 O 4- OH 

18 

O 4* H 2 O 2 ^ HOz 4" OH 

19 

OH + OH + M — H 2 O 2 4- M 

20 

OH 4- OH Hz 4- O 2 


Species : H 2 , O 2 , H 2 O, OH, H , O, HO 7 , H 2 O 2 and 
N 2 (inert) 

M is a third body (all species included) 


Table 2. Conditions for coaxial jet experiment 



Hz 

Air 

Mach No. 

1.0 

1.97 

Temperature 

545 K 

1180 K 

Pressure 

0.112 MPa 

0.107 MPa 

Velocity 

1774 m/s 

1380 m/s 

Of, 

1.0 

0.0 


0.0 

0.246 


0.0 

0.545 

f HiO 

0.0 

0.209 


Fuel injector diameter=0.00236 m 

Lip thickness=0. 000725 m 

Nozzle diameter(air flow) =0.0 1778 m 
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Table 3. Conditions for Burrows- Kurkov 
Experiment 



H 2 

Air 

Mach No. 

1.0 

2.44 

Temperature 

254 K 

1270 K 

Pressure 

0.1 MPa 

0.1 MPa 

Velocity 

1216 m/s 

1764 m/s 


1.0 

0.0 

fo. 

0.0 

0.258 

fa? 

0.0 

0.486 

t/fjO 

0.0 

0.256 


Fuel injector height=0.004 m 
Duct height at inlet=0.0938 m 
Duct height at exit=0.1048 m 



Uajj. = 1380m/s T air =1180K = 107000 N/m 2 

= 1774 m/s Ttt = 545 K P TT = 112000 N/m 2 

2 . H 2 H 2 

Fig. 1 Coaxial jet case : schematic 
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a) Total temperature 


b) Mach number 


30 


Fig. 6 Burrows-Kurkov ExpL * Temperature and Mach number at exit 



P ref = 91700 N/m 2 

Fig. 7 Burrows-Kurkov expt. : wall pressure 
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